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Accomplishments 

The  work  carried  out  under  this  contract  may  be  subdivided  according  to  the  following 
topics;  1)  Development  of  Implicit  ALE  Navier-Stokes  Solvers; 

2)  Implementation  of  Turbulence  Models; 

3-)  Development  of  Gridding  Techniques  for  3-D  Viscous  Flows; 

4)  Demonstration  Calculation  and  Results. 

The  description  of  the  main  accomplishments  of  the  present  work  are  listed  according 
to  these  topics  in  the  following. 

1.  IMPLICIT  ALE  NAVIER-STOKES  SOLVERS 

The  work  performed  in  the  area  of  flow  solvers  may  be  subdivided  into  two  parts; 

a)  Improvement  of  Edge-Based  Solvers; 

b)  Implicit  Solvers  for  the  Navier-Stokes  eqiiations. 

1.1  Edye-Baaed  Solvers:  The  switch  from  element-based  to  edge-based  solvers 
allowed  a  number  of  improvements  in  the  performance  and  accxiracy  of  the  flow- 
codes.  When  developing  these  new  codes,  we  incorporated  all  the  coding  lessons  that 
we  learned  over  the  years.  The  latest  code,  FEFL096,  nms  at  a  sustained  rate  of 
115  MFlops  on  the  CRAY-YMP,  and  has  sipiificantly  less  Flops  per  update  than  the 
previous  code  (FEFL052).  At  the  same  time,  indirect  addressing  (i/a)  costs  were  re¬ 
duced  by  a  factor  of  7.3.  This  reduction  in  i/a  was  achieved  by  using  the  edge-based 
data  structure  (1.57),  computing  the  fluxes  ‘on  the  fly’  from  the  unknowns  (2.14)  and 
using  superedges  (2.18)  [FI].  The  code  without  superedges  runs  about  25%  faster  on 
the  CRAY-YMP  than  the  code  using  usual  edges.  We  also  excercised  FEFL096  in 
parallel  on  the  CRAY-C90,  and  simply  using  autotasking,  i.e.  CRAY-preprocessing, 
timed  a  remarkable  3.87  speed-up  on  4  processors  for  some  large  runs. 
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1.2  Implicit  Solvers:  The  use  of  implicit  solvers  was  recognized  as  being  essential 
for  the  f2ist  solution  of  the  compressible  Navier-Stokes  equations.  Therefore,  we  devel¬ 
oped  implicit  iterative  solvers.  The  solvers  are  based  on  a  linearized  Euler  backward 
time-advancement  of  the  Navier-Stokes  equations.  The  resulting  non-linear  system  of 
equations  is  then  solved  using  a  preconditioned  GMRES  solver.  Both  block-diagonal 
and  incomplete  LU  decompositions  were  tried.  While  incomplete  LU  performs  bet¬ 
ter  than  block-diagoned,  it  also  requires  a  lot  more  memory,  and  is  less  amenable  to 
massive  parallel  architectures.  For  this  reason,  we  are  currently  investigating  other 
preconditioners.  More  details  of  these  techniques  are  given  in  [4],  which  is  included  in 
the  current  report  as  Appendix  1. 

2.  TURBULENCE  MODELS 

The  algebraic  Baldwin-Lomax  turbulence  model  [Tl]  was  implemented  in  the  3-D  flow 
solver.  This  is  the  simplest  way  to  add  the  effects  of  turbulence.  The  model  is  straight¬ 
forward  to  implement  as  the  only  difference  in  tho  flow  equations  is  the  increase  of 
viscosity.  The  model  represents  the  inner  and  outer  pairts  of  a  boundary  layer  zis  fol¬ 
lows: 

Inner  part: 

fit  =  p\uj\  [iy(l  -  exp{-y/A))\^ 

Outer  part: 

Pt  =  OtC cppYrnax^ mazy  > 

with 

^(y)  =  ykl(i  -  €xp(-y/A))  . 

The  exponential  factor  in  the  inner  part  is  the  Van  Driest  damping  factor  which  matches 
the  damping  of  the  wall.  In  the  second  equation,  Fmax  is  computed  along  a  normal 
to  the  wall.  This  is  pzurticularly  easy  when  dealing  with  structured  grids,  but  requires 
some  effort  for  unstructured  grids.  The  required  data  structures  were  discussed  by 
Rostand  [T2]  for  2-D.  As  far  as  we  can  see,  this  is  the  first  time  a  fully  unstructured 
3-D  Baldwin-Lomax  implementation  has  been  attempted.  One  complete  evaluation  of 
the  turbulent  viscosity  requires  the  following  transfer  of  information: 
a)  Vorticity  from  elements  or  points  in  the  mesh  to  the  appropriate  normals  to  the  walls. 
Along  each  normal,  the  vorticity  |u;|  is  required  to  evaluate  the  turbulent  viscosity.  This 
transfer  of  vorticity  is  accomplished  with  a  linked  list  of  intersections  of  normals  to  the 
wall  with  elements  of  the  mesh.  This  list  is  constructed  by  starting  from  the  surface  and 
moving  along  the  normal.  All  that  is  required  to  do  so  is  a  list  of  elements  surrounding 
elements. 
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b)  Transfer  of  the  turbulent  viscosity  from  the  normals  to  the  walls  to  the  elements  or 
points  of  the  mesh.  In  the  present  case,  we  transfer  to  points.  This  not  only  reduces  the 
transfers  required  (there  are  less  points  than  elements  in  a  mesh),  but  also  introduces  a 
benefici2il  smoothing  effect  at  element  level.  For  each  point  in  the  mesh  close  to  wetted 
surfaces,  we  find  the  closest  normals  surrounding  it,  and  then  the  two  closest  points 
along  each  of  the  normals.  Thus,  a  point  in  the  mesh  assembles  the  turbulent  viscosity 
from  four  points  along  the  normals  of  a  wetted  surface.  The  search  for  the  closest 
points  is  done  using  quad- trees  [T3].  In  the  case  of  the  mixing  of  several  boundary 
layers,  /it  is  computed  from  the  contribution  of  each  wall  weighted  by  its  distance  to 
the  point. 


3.  GRIDDING  FOR  3-D  VISCOUS  FLOWS 

The  difficulty  of  gi'idding  complex  geometries  for  the  simulation  of  flows  using  the 
Navier-Stokes  equations  -  i.e.  including  the  effects  of  viscosity  and  the  associated 
boundary  or  mixing  layers  -  increases  not  only  with  the  geometric  complexity  of  the  do¬ 
main  to  be  gridded,  but  also  with  the  Reynolds-number  of  the  flow.  For  high  Reynolds- 
numbers,  the  proper  discretization  of  the  very  thin,  yet  important  boundeiry  or  mixing 
layers  requires  elements  with  aspect  ratios  well  in  excess  of  1:1000.  This  requirement 
presents  formidable  difficulties  to  general,  ‘black-box’  unstructured  grid  generators. 
These  difficulties  can  be  grouped  into  two  main  categories: 

a)  Amount  of  Manual  Input:  In  most  unstructured  grid  generators,  the  desired  spatial 
distribution  of  element  size  and  shape  is  given  by  some  form  of  background  grid  or 
sources  [Gl,G2].  This  seems  natural  within  an  adaptive  context,  as  a  given  grid, 
combined  witn  a  suitable  error  indicator /estimator,  can  then  be  used  as  a  background 
grid  to  generate  an  even  better  grid  for  the  problem  at  hand.  Consider  now  trying 
to  generate  from  mzmual  input  a  first  grid  that  achieves  stretching  ratios  in  excess  of 
1:1,000.  The  amount  of  background  gridpoints  or  sources  required  will  be  proportional 
to  the  curvature  of  the  objects  immersed  in  the  flowfield.  This  implies  an  enormous 
amount  of  manual  labor  for  general  geometries,  rendering  this  approach  impractical. 

b)  Loss  of  Control:  Most  unstructured  grid  generators  introduce  a  point  or  element  at  a 
time,  checking  the  surrounding  neighbourhood  for  compatibility.  These  checks  involve 
Jacobians  of  elements  and  their  inverses,  distance-functions,  and  other  geometrical 
operations  that  involve  multiple  products  of  coordinate  differences.  It  is  not  difficult  to 
see  that  as  the  stretching-ratio  increases,  round-off  errors  can  become  a  problem.  For 
a  domain  spanning  1,000m  (mesh  axound  a  C-17),  with  a  minimum  element  length  at 
the  wing  of  0.0001m  across  the  boundary  layer  and  0.05m  along  the  boundary  layer, 
and  a  maximum  element  length  of  20m  in  the  farfield,  the  ratio  of  element  volumes  is 
of  the  order  of  3  *  10~^^.  Althoughthis  is  well  within  reach  of  the  10“^® -accuracy  of 
64-bit  arithmetic,  element  distortion  and  surface  singularities,  as  well  as  loss  of  control 
of  element  shape  can  quickly  push  this  ratio  to  the  limit. 

A  number  of  semi-automatic  grid  generators  have  been  devised  in  the  past.  The  most 
common  way  to  generate  meshes  suitable  for  Navier-Stokes  calculations  for  complex 
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geometries  is  to  employ  a  structured  or  semi-structuredmesh  close  to  wetted  surfaces 
or  wakes  [G3-G10].  This  ‘Navier-Stokes’  region  mesh  is  then  linked  to  an  outer  un¬ 
structured  grid  that  covers  the  ‘inviscid’  regions.  In  this  way,  the  geometric  complexity 
is  solved  using  unstructured  grids  and  the  physiczJ  complexity  of  near-wall  or  wake 
regions  is  solved  by  semi-structured  grids.  This  approach  has  proven  very  powerful  in 
the  past,  as  evidenced  by  many  examples. 

A  recurring  problem  in  all  of  these  approaches  has  been  how  to  link  the  semi-structured 
mesh  region  with  the  unstructured  mesh  region.  We  developed  a  new,  general  technique 
to  solve  this  problem.  The  design  criteria  for  the  new  grid  generation  strategy  may  be 
summarized  as  follows: 

-  The  geometric  flexibility  of  the  unstructured  grid  generator  should  not  be  com¬ 
promised  for  Navier-Stokes  meshes.  This  implies  using  unstructured  grids  for  the 
surface  discretization. 

-  The  memual  input  required  for  a  desired  Navier-Stokes  mesh  should  be  as  low  as 
that  used  for  the  Euler  czise.  In  the  present  C2ise,  this  requirement  is  solved  by 
specifying  at  the  points  of  the  background  grid  the  boundary  layer  thickness  and 
the  geometric  progression  normal  to  the  surface. 

-  The  generation  of  the  semi-structured  grid  should  be  fast.  Experience  shows  that 
usually  more  than  half  of  the  elements  of  a  typical  Navier-Stokes  mesh  cire  located 
in  the  boundary-layer  regions.  This  requirement  is  met  by  constructing  the  semi- 
structured  grids  with  the  same  normals  as  encountered  on  the  surface,  i.e.  without 
recurring  to  smoothing  procedures  as  the  semi-structured  mesh  is  advanced  into 
the  field  [G9,Gll]. 

-  The  element  size  and  shape  should  vary  smoothly  when  going  from  the  semi- 
structured  to  the  fully  unstructured  mesh  regions.  How  to  accomplish  this  is  the 
main  topic  of  this  paper,  and  is  deteuled  in  subsequent  sections. 

-  The  grid  generation  procediure  should  avoid  all  of  the  problems  typically  associated 
with  the  generation  of  Navier-Stokes  meshes  for  regions  with  high  surface  curva¬ 
ture:  negative  or  deformed  elements  due  to  converging  normals,  and  elements  that 
get  too  large  due  to  diverging  normals  at  the  surface.  In  order  to  circumvent  these 
problems,  the  same  techniques  which  are  used  to  achieve  a  smooth  matching  of 
semi-structured  emd  unstructured  mesh  regions  are  used. 

Given  these  design  criteria,  as  well  cis  the  approaches  used  to  meet  them,  the  present 
grid  generation  algorithm  cam  be  summarized  as  follows: 

M.l  Given  a  surface  definition  and  a  background  grid,  generate  a  surface  triauigulation 
using  an  unstructured  grid  generator. 

M.2  FVom  the  surface  triangulation,  obtain  the  surface  normeds. 

M.3  Smooth  the  surface  normals  in  several  passes  in  order  to  obtain  a  more  uniform 
mesh  in  regions  with  high  surface  curvature. 

M.4  Construct  a  semi-structured  grid  with  the  information  provided  by  the  background 
grid  and  the  smoothed  normals. 
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M.5  Examine  each  element  in  this  semi-structured  region  for  size  and  shape;  remove 
all  elements  that  do  not  meet  certain  specified  quality  criteria. 

M.6  Examine  whether  elements  in  this  semi-structured  region  cross  each  other;  if  so, 
keep  the  smaller  elements  and  remove  the  larger  ones,  until  no  crossing  occurs. 
M.6  Ex8unine  whether  elements  in  this  semi-structured  region  cross  boundaries;  if  so, 
remove  the  crossing  elements. 

M.7  Mesh  the  as  yet  ‘empty’  regions  of  the  computational  domain  using  the  background 
grid  and  the  unstructured  grid  generator. 

The  main  areas  of  work  were: 

a)  Smoothing  of  Surface  Normals:  This  implies  obtaining  the  smoothed  normals 
quickly  (i.e.  less  than  20  passes  over  the  surface  mesh),  and  defining  proper  bound¬ 
ary  conditions  for  the  normals  in  order  to  avoid  problems  at  ridges,  intersections, 
etc. 

b)  Prism  Generation:  When  generating  tetrahedra  from  prisms,  certeiin  compatibility 
criteria  must  be  met  [G14].  We  developed  a  very  fast  compatibility  algorithm  that 
was  found  to  converge  in  less  than  3  passes  over  the  surface  mesh. 

c)  Fast  Proximity  Finders  for  Crossed  Elements/Points:  In  order  to  speed  up  the 
checking  of  bad  elements  (negative,  crossing,  etc.),  a  series  of  data  structures  had 
to  be  developed  and  implemented.  The  main  data  structures  used  are  Octrees  and 
Linked  Lists.  The  main  filtering  strategies  to  reduce  the  work  even  further  are 
cones  of  visibility  and  distance  functions. 

d)  Initial  Front  for  Unstructured  Grid  Region:  Once  the  elements  have  been  marked 
for  deletion,  the  initisd  front  for  the  imstructured  grid  generation  of  the  remaining 
‘empty’  region  of  space  has  to  be  obteuned.  The  main  work  here  was  to  obtain 
this  initial  front  taking  into  consideration  the  boundary  arrays  required  by  the 
flow-solver  later  on,  i.e.  to  minimize  user-intervention  as  much  as  possible. 

More  details  of  the  algorithm,  as  well  as  several  example  grids  computed  with  it,  are 
given  in  [3],  which  is  included  in  the  report  as  Appendix  2. 

4.  DEMONSTRATION  CALCULATIONS  AND  RESULTS 

4.1  Unstructured  Grid/Remeshine  Transient  3-D  Runs:  With  the  developed  tools  we 
performed  several  2-D  and  3-D  test  runs.  The  steady  2-D  and  3-D  results  are  summa¬ 
rized  in  [1],  which  is  included  here  as  Appendix  3.  Pilot/Seat  ejection  from  an  F-16 
fighter  at  supersonic  speeds  was  shown  in  [2],  which  is  included  here  as  Appendix  4. 

5.  PUBLICATIONS 

All  of  the  developments  listed  above  were  reported  extensively  in  the  literature.  The 
main  papers  published  are  listed  in  chronological  order: 
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[1]  H.  Luo,  J.D.  Baum,  R.  Lohner  and  J.  Cabello  -  Adaptive  Edge-Based  Finite  Ele¬ 
ment  Schemes  for  the  Euler  and  Navier-Stokes  Equations;  AL4A-93-0336  (1993). 

[2]  J.D.  Baum  and  R.  Lohner  -  Numerical  Simulation  of  Pilot/Seat  Ejection  from  an 
F-16;  AIAA-93-0783  (1993). 

[3]  R.  Lohner  -  Matching  Semi-Structured  and  Unstructured  Grids  for  Navier-Stokes 
Calculations;  AIAA-93-3348-CP  (1993). 

[4]  H.  Luo,  J.D.  Baum,  R.  Lohner  and  J.  Cabello  -  An  Imphcit  Three-Dimensionsd 
Finite  Element  Solver  for  Unstructured  Meshes;  pp.  1027,1028  in  Proc.  11th  AIAA 
CFD  Conf.  ,  Orlando,  FL,  July  (1993). 

[5]  H.  Luo,  J.D.  Baum,  R.  Lohner-  Numerical  Solution  of  the  Euler  Equations 
for  Complex  Aerodynamic  Configurations  Using  An  Edge-Based  Finite  Element 
Scheme;  AIAA-93-2933  (1993). 

6.  CONCLUSIONS  AND  OUTLOOK 

We  have  made  major  steps  towards  the  development  of  a  CFD  capability  to  compute 
viscous  3-D  compressible  fiows  with  moving  bodies. 

In  the  future,  we  plan  to  expzmd  the  developed  capabilities  as  follows: 

a)  Development  of  Flow  Solvers  for  high  Re-number  Viscous  Flows:  We  will  continue 
to  develop  our  implicit  Navier-Stokes  solvers.  In  particular,  we  plan  to  extend  the 
linelet-concept  to  the  fully  coupled  linear  equation  systems  that  arise  for  implicit 
Navier-Stokes  solvers. 

b)  Turbulence  Models:  We  will  incorporate  the  k  —  c  model  into  the  implicit  Navier- 
Stokes  solver. 

c)  Development  of  suitable  gridding  algorithms  for  high  Re-number  viscous  fiows: 
We  plan  to  improve  our  present  capability  for  automatically  gridding  problems 
involving  high  Re-number  viscous  flows. 

d)  Development  of  suitable  error  indicators  for  high  Re-number  viscous  flows:  The 
idea  here  is  to  develop  error  indicators  that  sense  were  to  refine  boundary  layers 
and  shear  layers,  and  that  work  even  for  highly  stretched  grids. 

e)  Improvements  in  movie-making  capabilities:  The  main  aim  of  this  improvement  is 
to  be  able  to  run  the  movie  on  the  workstation  before  going  to  the  movie-making 
center.  At  the  same  time,  we  must  be  able  to  compress  the  information  to  an 
extent  that  a  3min  movie  can  be  stored  effortlessly  on  a  small  disk.  This  will  be 
accomplished  through  image  compression  algorithms. 

f)  Further  test  runs:  we  plan  to  look  for  available  experimentzJ  and/or  numerical 
data  in  the  literature  in  order  to  set  up  runs  to  test  the  algorithms  developed. 
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INTRODUCTION 

Significant  progress  has  been  made  in  recent 
years  in  developing  numerical  algorithms  for  the  so¬ 
lution  of  the  compressible  Euler  and  Navier-Stokes 
equations  on  unstructured  grids.  Most  efforts  have 
been  focused  on  improving  the  spatial  discretization 
operator  which  has  reached  a  high  degree  of  sophisti¬ 
cation.  Usually,  explicit  time  integration,  such  as  the 
multi-stage  Runge-Kutta  scheme  has  been  used  to  get 
a  steady  state  solution.  In  general,  explicit  schemes 
sue  easy  to  implement  and  vectorize  and  require  only 
limited  memory  storage.  However,  for  large-scale 
problems  and  especially  for  solution  of  the  Navier- 
Stokes  equations,  the  rate  of  convergence  slows  dra¬ 
matically.  To  speed  up  the  convergence  rate,  an  im¬ 
plicit  temporal  discretization  is  required. 

The  objective  of  this  research  is  to  develop  an  im¬ 
plicit  3D  finite  element  algorithm  for  the  solution  of 
the  compressible  Euler  and  Navier-Stokes  equations 
on  unstructured  meshes.  Numerical  results  for  both 
inviscid  and  viscous  flows  are  presented  to  demon¬ 
strate  the  performance  of  the  proposed  scheme. 

NUMERICAL  SCHEME 

The  governing  equations  are  integrated  in  time 
using  an  Euler  implicit  differencing  scheme.  The  re¬ 
sulting  large  nonsymmetric  linear  system  of  equations 
is  solved  by  using  the  preconditioned  GMRES  algo¬ 
rithm.  The  preconditioner  used  in  this  work  is  a  block 
.jomplete  LU  factorization  with  zero  fill-in.  Left, 
right,  and  symmetric  preconditioners  are  investigated 
and  examined.  Spatial  discretization  is  achieved  by 
using  an  edge-based  finite  element  scheme  [1].  Roe’s 
flux-difference  splitting  is  used  for  spatial  discretiza¬ 
tion  of  the  inviscid  flux  terms.  A  MUSCL  approach 
is  used  to  achieve  higher-order  accuracy.  The  viscous 
flux  terms  are  evaluated  using  second  order  accurate 
central  differences. 

Results  obtained  during  this  investigation  indi¬ 
cate  that  the  treatment  of  boundary  conditions  is 
extremely  important  to  the  success  of  an  implicit 
scheme.  When  boundary  conditions  are  treated  ex¬ 
plicitly,  only  a  very  li  nited  CFL  number  can  be  used. 
In  order  for  the  implicit  scheme  to  be  stable  for 
high  CFL  numbers,  the  boundary  condition  must  be 
treated  implicitly.  In  the  present  work.  Roe’s  approx¬ 
imate  Riemann  solver  was  also  extended  to  treat  the 


boundary  points  similarly  to  the  interior  points.  This 
procedure  provides  a  boundary  point  treatment  that 
is  completely  compatible  and  consistent  with  the  inte¬ 
rior  point  differencing  scheme.  The  boundary  condi¬ 
tions  are  then  linearized  consistently,  and  are  included 
in  the  left-hand-side  coefficient  matrix. 

NUMERICAL  RESULTS 
Due  to  space  limitation,  results  are  presented 
only  for  an  inviscid  transonic  flow  and  a  laminar  vis¬ 
cous  flow.  The  first  test  case  is  the  well  known  Ni’s 
test  case.  It  is  an  inviscid  flow  in  a  channel  with  a 
10%  thick  circular  bump  on  the  bottom.  Inlet  Mach 
number  is  0.675.  This  is  a  3D  simulation  of  a  2D 
flow.  The  mesh,  which  contains  13,891  grid  points. 
68,097  element  and  4,442  boundary  points,  is  depicted 
in  Fig. la.  Fig.  lb  displays  the  computed  pressure  con¬ 
tours  in  the  flow  field  The  Mach  number  distribution 
on  lower  wall  is  shown  in  Fig.lc.  Fig. Id  displays  a 
comparison  of  convergence  histories  between  explicit 
scheme  and  implicit  scheme  with  left,  right,  and  sym¬ 
metric  ILU  preconditioner,  respectively.  The  explicit 
scheme  results  were  obtsuned  using  three  stage  Runge- 
Kutta  scheme  with  implicit  residual  smoothing  and  a 
CFL  number  of  4.  The  implicit  scheme  results  were 
obtained  using  a  CFL  number  of  100,000  with  a  max¬ 
imum  number  of  Krylov  space  of  10  without  restart. 
Contrary  to  the  results  obtaind  in  [2],  the  left,  right 
and  symmetric  preconditioners  perform  equally  well. 

The  second  test  case  involves  a  3D  simulation  of 
a  2D  laminar  flow  past  a  flat  plate  at  a  Mach  num¬ 
ber  of  0.5  and  a  chord  Reynolds  number  of  10,000. 
The  mesh,  shown  in  Fig.2a.  contains  81,885  elements, 
15,694  points,  and  3,774  boundary  points.  The  com¬ 
puted  Mach  number  contours  and  velocity  vectors  in 
the  flow  field  are  depicted  in  Fig. 2b  and  2c.  Fig. 2d 
shows  the  comparison  of  the  Blasius  velocity  profile 
and  the  computed  velocity  profiles  as  scaled  by  the 
Blasius  similarity  law  at  the  exit. 
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ABSTRACT 

A  new  gridding  technique  for  Navier-Stokes  calculations  in¬ 
volving  complex  geometries  is  presented.  This  technique 
is  based  on  a  combination  of  semi-structured  and  unstruc¬ 
tured  meshing  techniques  that  accommodates  the  strengths 
of  these  two  approaches  while  avoiding  their  respective  weak¬ 
nesses.  The  technique  has  the  advantage  of  being  generally 
applicable,  yielding  one  single  unstructured  final  mesh  for  a 
given  computational  domain.  At  the  same  time,  the  prob¬ 
lems  usually  encountered  when  meshing  surfaces  with  high 
curvature  for  Navier-Stokes  calculations  are  avoided  auto¬ 
matically. 

1.  INTRODUCTION 

The  task  of  gridding  complex  geometries  for  the  simulation 
of  flows  using  the  Navier-Stokes  equations  -  i.e.  including 
the  effects  of  viscosity  and  the  associated  boundary  or  mix¬ 
ing  layers  •  is  encountered  conunonly  in  engineering  prac¬ 
tice.  The  difficulty  of  this  task  increases  not  only  with  the 
geometric  complexity  of  the  domain  to  be  gridded,  but  also 
with  the  Reynolds-number  of  the  flow.  For  high  Reynolds- 
numbers,  the  proper  discretization  of  the  very  thin,  yet  im¬ 
portant  boundary  or  mixing  layers  requires  elements  with 
aspect  ratios  well  in  excess  of  1:1,000.  This  requirement 
presents  formidable  difficulties  to  general,  ‘black-box’  un¬ 
structured  grid  generators.  These  difficulties  can  be  grouped 
into  two  main  categories: 

a)  Amount  of  Manual  Input:  In  most  unstructured  grid 
generators,  the  desired  spatial  distribution  of  element  size 
and  shape  is  given  by  some  form  of  background  grid  or 
sources  [1,2].  This  seems  natural  within  an  adaptive  con¬ 
text,  as  a  given  grid,  combined  with  a  suitable  error  in¬ 
dicator/estimator,  can  then  be  used  as  a  background  grid 
to  generate  an  even  better  grid  for  the  problem  at  hand. 
Consider  now  trying  to  generate  from  manual  input  a  first 
grid  that  achieves  stretching  ratios  in  excess  of  1:1,000.  The 
amount  of  background  gridpoints  or  sources  required  will  be 
proportional  to  the  curvature  of  the  objects  immersed  in  the 
flowfield.  This  implies  an  enormous  amount  of  muual  labor 
for  general  geometries,  rendering  this  appro^udl  impractical. 


ement  length  of  20m  in  the  farfield,  the  ratio  of  element 
volumes  is  of  the  order  of  3  *  10“'^.  Although  this  is  well 
within  reach  of  the  10“^*-accuracy  of  64-bit  arithmetic,  el¬ 
ement  distortion  and  surface  singularities,  as  well  as  loss  of 
control  of  element  shape  can  quickly  push  this  ratio  to  the 
limit. 

Given  these  difficulties,  it  is  not  surprising  that  at  present, 
there  does  not  exist  a  ‘black-box’  unstructured  (or  struc¬ 
tured,  for  that  matter)  grid  generator  that  can  produce  ac¬ 
ceptable  meshes  with  such  high  aspect  ratio  elements.  The 
demand  for  Navier-Stokes  calculations  in  or  past  complex  ge¬ 
ometries  being  great,  a  number  of  semi-automatic  grid  gener¬ 
ators  have  been  devised.  The  most  common  way  to  generate 
meshes  suitable  for  Navier-Stokes  calculations  for  complex 
geometries  is  to  employ  a  structured  or  semi-structured  mesh 
close  to  wetted  surfaces  or  wakes  [3-llj.  This  ‘Navier-Stokes’ 
region  mesh  is  then  linked  to  an  outer  unstructured  grid  that 
covers  the  ‘inviscid’  regions.  In  this  way,  the  geometric  com¬ 
plexity  is  solved  using  unstructured  grids  and  the  physical 
complexity  of  near-wall  or  wake  regions  is  solved  by  semi- 
stnictured  grids.  This  approach  has  proven  very  powerful  in 
the  past,  as  evidenced  by  many  examples. 

The  meshes  in  the  semi-structured  region  can  be  constructed 
to  be  either  quads/bricks  [3-9]  or  triangles/prisms  [10,11]. 
The  prisms  can  then  be  subdivided  into  tetrahedra  if  so 
desired.  For  the  inviscid  (unstructured)  regions,  multi¬ 
block  approaches  have  been  used  for  quads/bricks  [5-8], 
and  advancing  front  [13-15],  Voronoi  [16,17]  and  modified 
quadtree/octree  apprortches  [18]  for  triangles/tetrrdiedra. 

A  recurring  problem  in  all  of  these  approaches  has  been  how 
to  link  the  semi-structured  mesh  region  with  the  unstruc¬ 
tured  mesh  region.  Some  solutions  put  forward  have  consid¬ 
ered: 

-  Overlapped  Structured  Grids:  these  are  the  so-called 
chimera  grids  [8,9],  that  have  become  popular  for  Navier- 
Stokes  calculations  of  complex  geometries;  as  the  grid- 
points  of  the  various  meshes  do  not  coincide,  they  allow 
great  flexibility  and  are  easy  to  construct,  but  the  so¬ 
lution  has  to  be  interpolated  between  grids,  which  ma\ 


h)  Loss  of  Control;  Most  unstructured  grid  generators  intro¬ 
duce  a  point  or  element  at  a  time,  checking  the  surrounding 
neighbourhood  for  compatibility.  These  checks  involve  Jaco- 
bians  of  elements  Md  their  inverses,  distance-functions,  and 
other  geometrical  operations  that  involve  multiple  products 
of  coordinate  differences.  It  is  not  difficult  to  see  that  as 
the  stretching-ratio  increases,  round-off  errors  can  become 
a  problem.  For  a  domain  spanning  1,000m  (mesh  around  a 
Boeing-747),  with  a  minimum  element  length  at  the  wing  of 
less  than  0.01mm  across  the  boundary  layer  and  0.05m  along 
^^e^hotmdary  layer  and  along  the  wing,  and  a  maximum  el- 


lead  to  higher  CPU  cost  and  a  deterioration  in  solutior 
quality. 

Overlapped  Structured /Unstructured  Grids:  in  this  cast 
the  overlap  zone  can  be  restricted  to  one  cell,  with  th< 
points  coinciding  exactly,  so  that  there  are  no  interpola 
tion  problems  [3,4]. 

Delaunay  Triangulation  of  Points  Generated  by  Alge 
braic  Grids:  in  this  case  several  structured  grids  are  gen 
erated.  and  their  spatial  mapping  functions  are  stored 
the  resulting  cloud  of  points  is  then  gridded  using  De 


^^Pynght  1993  by  the  author.  Published  by  the  AIAA  with  permission. 
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Uunay  triangulation  techniques  [16]. 

I 

Although  some  practical  problems  have  been  solved  by  these 
approaches,  they  can  not  be  considered  general,  as  they  suf¬ 
fer  from  the  following  constreunts: 

-  The  first  two  approaches  require  a  very  close  link  be¬ 
tween  solver,  grid  generator  and  interpolation  techniques 
to  achieve  good  results;  from  the  standpoint  of  general¬ 
ity,  such  a  close  link  between  solver,  grid  generator  and 
interpolation  modules  is  undesirable. 

-  Another  problem  associated  with  the  first  two  ap¬ 
proaches  is  that  at  concave  comers,  negative  (i.e.  folded) 
or  badly  shaped  elements  may  be  generated.  The  usual 
recourse  is  to  smooth  the  mesh  repeatedly,  or  use  some 
other  device  to  introduce  ellipticity  [10,12].  These  ap¬ 
proaches  tend  to  be  CPU  intensive,  and  require  consid¬ 
erable  expertise  from  the  user.  Therefore,  they  csin  not 
be  considered  general  approaches. 

-  The  third  case  requires  a  librtuy  of  algebraic  grids  to 
mesh  individual  cases,  and  can  therefore  not  be  consid¬ 
ered  a  general  tool.  However,  it  has  been  used  exten¬ 
sively  for  important  specialized  applications,  e.g.  single 
or  multi-element  airfoil  flows  [16]. 

The  present  research  effort  is  directed  towards  generality  and 
ease  of  software  expandability  and  maintainability.  There¬ 
fore,  we  strive  to  generate  a  single  unstructured  mesh  consist¬ 
ing  of  triangles/tetrahedra.  This  mesh  can  then  be  consid¬ 
ered  completely  independent  of  flow  solvers,  and  neither  re¬ 
quires  any  interpolation  or  other  transfer  operators  between 
grids,  nor  the  storage  of  mapping  functions. 

The  remainder  of  the  paper  is  divided  as  follows:  The  de¬ 
sign  criteria  used  and  the  new  grid  generation  strategy  are 
outlined  in  Section  2.  Removal  criteria  are  discussed  in  Sec¬ 
tion  3.  Section  4  uescribes  the  formation  of  tetrafaedra  from 
prismatic  elements.  Several  examples  of  gridded  configura¬ 
tions  are  shown  in  Section  5.  Finally,  some  conclusions  are 
drawn  in  Section  6,  tmd  further  extensions  are  considered. 

2.  DESIGN  CRITERIA  AND  ALGORITHM 

The  design  criteria  for  the  new  grid  generation  strategy  may 
be  summarized  as  follows: 

-  The  geometric  flexibility  of  the  unstructured  grid  genera¬ 
tor  should  not  be  compromised  for  Navier-Stokes  meshes. 
This  implies  using  unstructured  grids  for  the  surface  dis¬ 
cretization. 

-  The  manual  input  required  for  a  desired  Navier-Stokes 
mesh  should  be  as  low  as  that  used  for  the  Euler  case.  In 
the  present  case,  this  requirement  is  solved  by  specifying 
at  the  points  of  the  background  grid  the  boundary  layer 
thickness  and  the  geometric  progression  normal  to  the 
surface. 

-  The  generation  of  the  semi-structured  grid  should  be 
fast.  Experience  shows  that  usually  more  than  half  of 
the  elements  of  a  typical  Navier-Stokes  mesh  are  located 
in  the  boundary-layer  regions.  This  requirement  is  met 
by  constructing  the  semi-structured  grids  with  the  same 
normals  as  encountered  on  the  surface  (see  Figure  1), 


i.e.  without  recurring  to  smoothing  procedures  as  the 
semi-structured  mesh  is  advanced  into  the  field  [10,12]. 

-  The  element  size  and  shape  should  vary  smoothly  when 
going  from  the  semi-structured  to  the  fully  unstructured 
mesh  regions.  How  to  accomplish  this  is  the  main  topic 
of  this  paper,  and  is  detailed  in  subsequent  sections. 

-  The  grid  generation  procedure  should  avoid  all  of  the 
problems  typically  associated  with  the  generation  of 
Navier-Stokes  meshes  for  regions  with  high  surface  cur¬ 
vature:  negative  or  deformed  elements  due  to  converging 
normals,  and  elements  that  get  too  large  due  to  diverg¬ 
ing  normals  at  the  surface.  In  order  to  circumvent  these 
problems,  the  same  techniques  which  are  used  to  achieve 
a  smooth  matching  of  semi-structured  and  unstructured 
mesh  regions  are  used. 

Given  these  design  criteria,  as  well  as  the  approaches  used 
to  meet  them,  the  present  grid  generation  ailgorithm  can  be 
summarized  as  follows  (see  Figure  1): 

M.l  Given  a  surface  definition  and  a  background  grid,  gen¬ 
erate  a  surface  triangulation  using  an  unstructured  grid 
generator. 

M.2  From  the  surface  triangulation,  obtain  the  surface  nor¬ 
mals. 

M.3  Smooth  the  surface  normals  in  several  passes  in  order  to 
obtain  a  more  uniform  mesh  in  regions  with  high  surface 
curvature. 

M.4  Construct  a  semi-structured  grid  with  the  information 
provided  by  the  background  grid  and  the  smoothed  nor¬ 
mals. 

M.5  Examine  each  element  in  this  semi-structured  region  for 
size  and  shape;  remove  all  elements  that  do  not  meet 
certain  specified  quality  criteria. 

M.6  Examine  whether  elements  in  this  semi-structured  region 
cross  each  other;  if  so,  keep  the  smaller  elements  and 
remove  the  larger  ones,  until  no  crossing  occurs. 

.M.7  Examine  whether  elements  in  this  semi-structured  region 
cross  boundairies;  if  so,  remove  the  crossing  elements. 

M.8  Mesh  the  as  yet  ‘empty’  regions  of  the  computational 
domain  using  an  unstructured  grid  generator  in  combi¬ 
nation  with  the  desired  element  size  and  shape. 

3.  ELEMENT  REMOVAL  CRITERIA 

The  critical  element  of  the  matching  rdgorithm  described 
above  is  the  development  of  good  element  removal  crite¬ 
ria.  The  criteria  to  be  considered  sire:  element  size,  element 
shape,  element  overlap  and  element  crossing  of  boundary 
faces. 

3.1  Element  Size 

The  two  main  types  of  problems  encountered  in  semi- 
structured  grid  regions  that  are  related  to  element  size  are 
elements  that  are  either  too  large  or  negative  (folded).  These 
problems  originate  for  different  reasons,  and  will  therefore  be 
treated  separately. 
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3.1.1  Large  Elementa 

As  a  result  of  surface  normals  diverging  close  to  convex  sur¬ 
faces  very  large  elements  (as  compared  to  the  iisa-defined 
size  and  shape)  may  appear  in  the  semi-structured  mesh  re¬ 
gions.  The  situation  is  shown  diagrammatically  in  Figure  2. 
The  volume  of  each  element  in  the  semi-structured  mesh  re¬ 
gion  is  compared  to  the  element  volume  desired  by  the  user 
for  the  particular  location  in  space.  Any  element  with  a  vol¬ 
ume  greater  than  the  one  specified  by  the  user  is  marked 
for  deletion.  In  the  present  case,  the  desired  distribution 
of  element  size  and  shape  is  given  by  a  background  grid. 
Tree-search  algorithms  are  used  to  relate  the  information 
between  this  background  grid  and  a  particular  location  in 
space  (see  (13]  for  details). 

3.1.2  Negative  Elements 

As  a  result  of  folding  away  from  concave  surfaces,  elements 
with  negative  jacobians  may  appear.  The  situation  is  shown 
diagrammatically  in  Figure  3.  As  before,  the  element  vol¬ 
umes  are  computed.  All  elements  with  negative  volumes  are 
marked  for  deletion. 

We  have  observed  that  typically  the  elements  adjacent  to 
negative  elements  tend  to  be  highly  deformed.  Therefore,  we 
also  remove  all  elements  that  have  points  in  common  with 
negative  elements.  Obviously,  this  one-pass  procedure  can  be 
extended  to  several  passes,  i.e.  neighbours  of  neighbours,  etc. 
Our  experience  indicates,  however,  that  one  pass  is  sufficient 
for  most  cases. 

3.2  Element  Shane 

The  aim  of  a  semi-structured  mesh  close  to  a  wall  is  to  pro¬ 
vide  elements  with  very  small  size  normal  to  the  wail  and 
reasonable  size  along  the  wall.  Due  to  different  meshing  re¬ 
quirements  along  the  wall  (e.g.  comers,  separation  points, 
leading  and  trailing  edges  for  small  element  size,  other  re¬ 
gions  with  larger  element  size),  elements  that  are  longer  in 
the  direction  normal  to  the  wall  than  along  the  wall  may 
appear.  The  situation  is  shown  diagrsuiunatically  in  Fig¬ 
ure  4.  For  the  semi-structured  grids,  the  element  and  point 
numbering  can  be  assumed  as  known.  Therefore,  a  local  el¬ 
ement  analysis  can  be  performed  to  determine  whether  the 
element  has  side-ratios  that  are  consistent  with  boundary 
layer  gridding.  All  elements  that  do  not  satisfy  this  criterion 
are  marked  for  deletion. 

3.3  Overlapping  Elements 

Overlapping  elements  will  occur  in  regions  close  to  concave 
surfaces  with  high  curvature,  or  when  the  semi-structured 
grids  of  two  close  objects  overlap.  Another  possible  scenario 
is  the  overlap  of  the  semi-structured  grids  of  mixing  wakes. 
The  main  criterion  employed  is  to  keep  the  smaller  element 
whenever  an  overlap  occurs.  In  this  way,  the  small  elements 
close  to  surfaces  are  always  retained.  Straightforward  testing 
would  result  in  OflV,;)  operations  per  element,  where  N,i 
denotes  the  number  of  elements,  leading  to  a  total  number 
of  operations  of  0(fV*j).  By  using  quad/octrees  [13],  or  other 


suitable  data  structures  [19],  the  number  of  dements  tested 
can  be  reduced  significantly,  leading  to  a  total  number  of 
operations  of  0{N,i  log  Nti). 

For  quad/octrees,  the  complete  testing  procedure  would  look 
as  follows: 

-  Construct  a  quad/octree  for  the  points; 

-  Order  the  elements  according  to  decreasing  volume  (e.g. 
in  a  heap-list  [13]); 

-  Construct  a  linked  list  for  ail  the  elements  surrounding 
each  point; 

-  Loop  over  the  elements,  in  descending  volume,  testing: 

-  IF  the  element,  denoted  in  the  sequel  by  IBJBI,  has  not 
been  marked  for  deletion  bdore: 

-  Obtain  the  minimum/maximum  extent  of  the  coor¬ 
dinates  belonging  to  this  element; 

-  Find  from  the  quad/octree  all  points  falling  into  this 
search  region,  storing  them  in  a  list  LCLOPdtlCLOP); 

-  Find  all  the  unmarked  elements  with  smaller  vol¬ 
ume  than  lELBH  surrounding  the  points  stored  in 
LCLap(i:icLOF);  this  yields  a  list  of  close  elements 
LCU»(1:ICL0B>; 

-  Loop  over  the  elements  stored  in  LCLOl(l:flCLOE): 

-  IF  the  element  crosses  the  faces  of  mill 

or  is  inside  lELSX:  Mark  lELBi  for  deletion 

The  reason  for  looping  over  the  elements  according  to  de¬ 
scending  volumes  is  that  the  search  region  is  obtained  in  a 
natural  way,  i.e.  the  extent  of  the  element.  Looping  ac¬ 
cording  to  ascending  volumes  would  imply  guessing  search 
regions. 

As  negative  elements  could  lead  to  a  failure  of  this  test,  the 
overlap  test  is  performed  after  the  negative  elements  have 
been  identified  and  marked. 

3.4  Elements  Crossine  Boundary  Faces 

In  regions  where  the  distance  between  surfaces  is  very  small, 
the  crossing  of  boundary  faces  by  elements  from  the  semi- 
structured  region  is  likely  to  occur.  As  this  test  is  per¬ 
formed  after  the  element  crossing  tests  are  conducted,  the 
only  boundaries  that  need  to  be  treated  are  those  that  have 
no  semi-structured  grid  attached  to  it.  In  order  to  detect 
if  overlapping  occurs,  we  loop  over  the  surface  faces,  see¬ 
ing  if  any  element  crosses  it.  As  before,  straightforward 
testing  would  result  in  an  expensive  0{Nti  •  Nf)  procedure, 
where  Nf  denotes  the  number  of  boundary  faces.  By  us¬ 
ing  quad/octrees  [13],  this  complexity  can  be  reduced  to 
0{Nf  log  Nft)-  The  face-crossing  check  looks  essentially  the 
same  as  the  check  for  overlapping  elements,  and  its  explicit 
description  is  therefore  omitted. 

4.  SUBDIVISION  OF  PRISMS  INTO  TETRAHEDRA 

As  we  only  desire  to  work  with  tetrahedra,  the  prisms  formed 
by  extruding  the  surface  triangles  along  the  smoothed  nor¬ 
mals  must  be  subdivided.  This  subdivision  must  be  per¬ 
formed  in  such  a  way  that  the  diagonals  introduced  at  the 


reclaaguiar  of  the  prisms  match  across  prisms.  Given  gap  regions.  Application  of  the  removal  criteria  reduced  the 

that  a  priaa  caonot  be  subdivided  into  tetrahedra  in  any  number  of  elements  to  IELEll»60,022,  yielding  the  mesh  shown 

arbitrary  ww,  care  has  to  be  taken  when  choosing  these  di-  in  Figure  7c.  The  Anal  unstructured  mesh,  consisting  c, 

agonais.  Figure  5  illustrates  the  possible  diagonals  as  the  iELEM*St,es7  elements  is  shown  in  Figure  7d. 


base  ^ides  'x  the  prism  are  traversed.  One  can  see  that  in 
order  to  obtain  a  combination  of  diagonals  that  can  be  sub¬ 
divided  into  tetrahedra,  not  all  sides  of  the  triangular  base 
have  to  be  sp-down  or  down-up  as  one  traverses  the  sides. 
This  implies  that  the  sides  of  the  triangular  base  mesh  have 
to  be  marked  in  such  a  way  that  no  such  combination  oc¬ 
curs.  We  have  implemented  the  following  iterative  procedure 
to  arrive  at  valid  aide-combinations: 

D.O  Given; 

■  The  sides  of  the  surface  triangulation 

-  Tie  sides  of  each  surface  triangle 

-  'Hje  triangles  that  surround  each  surface  triangle 
D.l  DO;  li>op  over  the  surface  triangles 

-  IF  the  current  side-combination  is  not  valid: 

-  DO;  loop  over  the  sides  of  the  triangle 
-  IF  the  inversion  of  the  side/diagonal 
orientation  leads  to  an  allowed 
side-combination  in  the 
neighbour-triangle; 

-  Invert  the  side/diagonal  orientation 
•  Goto  next  element 
BIOIF 
EISDO 

ODIF 

EIDDO 

D.2  IF  any  unallowed  combinations  are  left:  GOTO  D.l 

This  procedure  works  well,  converging  in  at  most  two  passes 
over  the  surface  mesh  for  all  cases  tested  to  date. 

5.  EXAMPLES 

lilorked  Ghannel  with  Object:  The  surface  definition  is 
given  in  Figure  6a.  The  semi-structured  grid  generated  from 
the  surface  discretization,  consisting  of  »SLEii»7.7B4  element, 
is  shown  in  Figure  fib.  Observe  that  we  have  sharp  con¬ 
vex  and  concave  corners,  smoothed  normals,  negative  ele¬ 
ments  at  concave  corners,  and  element  overlap  between  semi- 
sltuciured  grids.  Application  of  the  removed  criteria  reduced 
the  number  of  elements  to  ■ELEH<4,638,  yielding  the  mesh 
shown  in  Figure  6c.  The  final  unstructured  mesh,  consisting 
of  IELEN>S,39T  elements  is  shown  in  Figure  6d,e.  Notice  the 
smooth  transition  between  the  semi-structured  and  unstruc¬ 
tured  grid  regions. 

Multi-F.lqment  Airfoil:  The  surface  definition  is  shown  in 
l•igure  7a.  The  semi-structured  grid  generated  from  the  sur¬ 
face  discretization,  is  shown  in  Figure  7b  and  consisted  of 
NELEM>E7,2S2  elements.  As  before,  the  normals  have  been 
smoothed,  but  deformed  elements  appear  due  to  surface  cur¬ 
vature.  Overlapped  elements  are  present  in  the  inter-airfoil 


rylinder  on  ^  Flat  Plate:  The  surface  definition  for  this  3- 
D  case  is  shown  in  Figure  8a.  The  surface  of  the  semi- 
structured  grid  generated  from  the  surface  triangtilation, 
which  consisted  of  iBLEMOsa.WO  elements,  is  shown  in  Fig¬ 
ure  8b.  Overlapping  and  negative  elements  are  clearly 
present  in  this  mesh.  Application  of  the  removal  criteria 
reduced  this  number  to  ISLB!«138.964  yieldmg  the  surface 
shown  in  Figure  8c.  The  removal  of  these  elements  required 
less  than  Imin  of  CPU  time  on  an  lBM-RISC-550  worksta<- 
tion.  The  surface  of  the  final  unstructured  mesh,  consisting 
of  IELSM<>213,979  elements  is  shown  in  Figure  8d.  The  CPU 
time  required  for  the  complete  mesh  was  less  than  lOmin 
on  an  IBM-RISC-550  workstation.  This  mesh  was  used  for 
an  incompressible  laminar  flow  simulation.  A  Blasius  pro¬ 
file  was  prescribed  at  the  entrance  plane,  and  the  Reynolds- 
number  based  on  the  cylinder  diameter  was  set  to  /le  =  100. 
Figures  8e,f  show  some  of  the  cross-sectional  meshes,  as  well 
as  the  solution  obtained. 

riOTeric  Missile:  The  surface  definition  for  this  3-D  case  is 
shown  in  Figure  9a.  The  surface  of  the  semi-stnictured  grid 
generated  from  the  surface  triangulation,  which  consisted 
of  reixii»94e,800  elements,  is  shown  in  Figure  9b.  Over¬ 
lapping  and  negative  elements  are  clearly  present  in  this 
mesh.  Application  of  the  removal  criteria  reduced  this  num¬ 
ber  to  IELEM»69T,638  yielding  the  surface  shown  in  Figure  9c. 
The  surface  of  the  final  unstructured  mesh,  consisting  of 
•m  r»»oo7  Ban  elements  is  shown  in  Figure  9d.  The  CPU 
time  required  for  the  complete  mesh  was  about  fifimin  on  an 
IBM-RISC-550  workstation. 

6.  CONCLUSIONS  AND  OUTLOOK 

A  new  gridding  technique  for  Navier-Stokes  calculations  in¬ 
volving  complex  geometries  has  been  presented.  The  tech¬ 
nique  is  based  on  a  combination  of  semi-structured  and 
unstructured  meshing  techniques  that  accommodates  the 
strengths  of  these  two  approaches  while  avoiding  their  re¬ 
spective  weaknesses.  The  technique  has  the  advantage  of 
being  generally  applicable,  yielding  one  single  unstructured 
final  mesh  for  a  given  computational  domain.  At  the  same 
time,  the  problems  usually  encountered  when  meshing  sur¬ 
faces  with  high  curvature  for  Navier-Stokes  calculations  we 
avoided  automatically.  The  technique  was  demonstrated  on 
several  examples  that  were  run  interactively  on  workstations, 
indicating  a  reasonable  speed  for  possible  use  within  an  adap¬ 
tive  remeshing  context.  Future  developments  will  center  on 
better  wake-gridding  capabilities,  automation  of  the  proce¬ 
dure  for  adaptive  remeshing,  and  the  extension  to  free  sur¬ 
faces  or  flexible  bodies  immersed  in  the  flowfield. 
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Figure  8;  Cylinder  on  a  Flat  Plate 
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Figure  9  Generic  Missile  Configuration 
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ABSTRACT 

This  paper  describes  recent  developments  of  high 
resolution  finite  element  schenrtes  for  the  solution  of 
the  unsteady  compressible  Euler  and  Navier-Stokes 
equations  on  unstructured  meshes.  These  finite  el¬ 
ement  algorithms  use  an  edge-based  data  structure, 
as  opposed  to  a  more  traditional  element-based  data 
structure.  The  advantage  of  using  such  an  edge-based 
data  structure  is  that  it  provides  a  unified  approach 
in  which  the  relation  between  centered  and  upwind 
schemes  becomes  apparent,  improves  the  efficiency  of 
the  algorithms,  and  reduces  the  storage  requirements. 

A  variety  of  numerical  schemes  using  such  edge-based 
data  structure,  ranging  from  Godunov  schemes  to 
centered  schemes  with  blended  dissipation,  is  pre¬ 
sented  and  discussed.  Adaptive  mesh  refinement  is 
then  added  to  these  solvers  to  enhance  the  solution 
accuracy  and  efficiency.  Various  numerical  results  for 
a  wide  range  of  flow  conditions,  from  subsonic  to  hy¬ 
personic  in  both  2D  and  3D,  are  presented  to  demon¬ 
strate  the  performance  and  versatility  of  the  proposed 
schemes. 

1.  INTRODUCTION 

In  recent  years,  significant  progress  has  been 
made  on  developing  numerical  algorithms  for  the  so¬ 
lution  of  the  compressible  Euler  and  Navier-Stokes 
equations.  The  use  of  unstructiued  meshes  for 
computational  fluid  dynamics  problems  has  become 
widespread  due  to  their  ability  to  discretize  arbitrar¬ 
ily  complex  geometries  and  the  ease  with  which  mesh 
adaption  can  be  carried  out  to  improve  the  solution. 
However,  any  numerical  schemes  based  on  unstruc¬ 
tured  meshes  require  a  storage  of  mesh  coimectivity 
information.  This  requirement  leads  to  an  increase 
of  computer  memory  and  the  use  of  indirect  address¬ 
ing  to  retrieve  nearest  neighbor  information,  which, 
in  turn,  implies  that  any  numerical  algorithms  will 
run  slower  on  unstructured  grids  than  on  structured 
grids.  To  reduce  indirect  addressing,  finite  element 
schemes  based  on  edge-based  data  structures  have 
been  introduced  [1-4].  In  addition,  more  sophisti¬ 
cated  data  structures  such  as  stars,  super  edges,  and 
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chains  have  recently  been  developed  [5].  The  use  of 
edge-based  data  structure  has  been  shown  to  yield  sig¬ 
nificant  computational  savings  for  three  diniensional 
problems. 

Extensive  research  has  been  performed  during 
the  last  few  years  on  upwind  algorithms  for  the  so¬ 
lution  of  the  Euler  and  Navier-Stokes  equations  on 
unstructured  meshes  [6-9].  A  significant  advantage  of 
any  upwind  discretization  is  that  it  is  naturally  dissi¬ 
pative,  as  compared  with  centraJ-difference  discretiza¬ 
tions,  and  consequently  does  not  require  any  problem- 
dependent  parameters  to  adjust.  ^  far,  all  the  up¬ 
wind  schemes  mq>lemented  as  either  node-centered  or 
cell-centered  discretizations  on  unstructured  meshes 
use  the  finite  volume  approach  and  the  control  vol¬ 
ume  must  be  constructed  first.  In  terms  of  computa¬ 
tional  efficiency,  node-centered  schemes  are  preferable 
to  their  cell-center  counterparts.  In  the  node-centered 
approach  [6,8],  the  control  volume  is  typically  taken 
to  be  part  of  the  neighboring  cells  that  have  a  ver¬ 
tex  at  that  node.  In  two  dimensions,  the  part  of  the 
cells  taken  is  determined  by  connecting  the  centroid 
of  the  cell  and  the  midpoints  of  the  two  edges  that 
share  the  node.  In  3-D,  the  part  of  the  cells  taken 
is  determined  by  a  surface  that  is  constructed  in  a 
similar  way,  a  somewhat  complicated  geometrical  pro¬ 
cess  in  three  dimensions.  The  switching  from  element 
to  edge-based  data  structure  enables  the  implemen¬ 
tation  of  upwind  schemes  trivial  and  straightforward 
in  the  context  of  finite  elements.  This  is  especially 
attractive  for  three  dimensional  problems,  as  there  is 
no  need  to  construct  control  volumes  explicitly  and 
geometrically. 

The  objective  of  this  paper  is  to  present  recently 
developed  high  accuracy  s^emes  on  uirstructured 
grids  using  an  edge-based  data  structure.  This  edge- 
based  data  structure  provides  a  unified  sq>proach  in 
which  the  link  between  centered  and  upwind  schemes 
becomes  apparent.  The  use  of  such  an  edge-based 
data  structure  not  only  improves  the  efficiency  of 
the  algorithms,  but  also  enables  a  straightforward 
implementation  of  upwind  schemes  in  the  context 
of  finite  element  methods.  A  variety  of  numerical 
schemes  using  the  edge-based  data  structure  is  pre- 


seated  and  the  performance  of  these  schemes  in  terms 
of  solution  accuracy  and  overall  computational  ef¬ 
ficiency  is  discussed.  Some  different  strategies  for 
the  discretization  of  the  viscous  terms  are  considered. 
An  approach  well  suited  for  use  with  an  edge-based 
data  structure  is  then  introduced  and  presented.  An 
H-refinement/coarsening  adaptive  scheme  is  imple¬ 
mented  in  these  schemes  to  enhance  the  solution  ac¬ 
curacy  and  efficiency.  Various  numerical  examples  for 
a  wide  range  of  flow  conditions,  from  subsonic  to  hy¬ 
personic  in  both  2D  and  3D,  are  presented  to  demon¬ 
strate  the  performance  and  versatility  of  the  proposed 
adgorithms. 

2.  GOVERNING  EQUATIONS 

The  Navier-Stokes  equations  governing  unsteady 
compressible  viscous  flows  can  be  expressed  in  the 
conservative  form  as 


ay  d& 

dt  dxj  dxj  ’ 


(2.1) 


where  the  summation  convention  has  been  employed 
and 


y  =  f  4,  y  P  f  puiu^i  pSij  \  , 

\pe/  \ujipe+p}J 


Here  p,p,e,Taad  k  denote  the  density,  pressure,  spe¬ 
cific  total  energy,  temperature  and  thermal  conduc¬ 
tivity  of  the  fluid,  respectively,  and  u,-  is  the  velocity 
of  the  flow  in  the  coordinate  direction  z,-.  This  set  of 
equations  is  completed  by  the  addition  of  the  equation 
of  state 


In  the  sequel,  we  assume  that  D  is  the  flow  do¬ 
main,  r  its  bounduy,  and  rij  the  unit  outward  normal 
vector  to  the  boundary.  The  following  bound  ^jry  con¬ 
ditions  have  to  be  added: 

On  the  solid  wall,  the  slip  condition  is  at  led  for 
inviscid  flow 

=  0  .  (2.6) 

For  viscous  flow,  no-slip  condition 

u,  =  0  ,  (2.7) 

and  either  isothermal  condition 

T=To  (2.8) 

where  To  is  the  total  temperature  or  adiabatic  condi¬ 
tion 

dT 

^n..=0  (2.9) 

could  be  imposed.  In  the  far-held,  a  characterbtic 

analysis  based  on  the  introduction  of  Riemann  invari¬ 
ants  for  one-dimensional  flow  normal  to  the  boundary 
is  used  to  determine  the  values  of  the  flow  variables. 
This  analysis  correctly  accounts  for  wave  propagation 
in  the  far  field,  which  is  important  for  rapid  conver¬ 
gence  to  steady-state  and  serves  as  a  non-reflecting 
boundary  condition  for  unsteady  application. 

3.  VARIATIONAL  FORMULATION  AND 
FINITE  ELEMENT  APPROXIMATION 

Let  T  be  a  trial  function  space  and  W  a  weighting 
function  space,  both  defined  to  consist  of  all  suitably 
snKx>tb  functions.  An  equivalent  variational  formula¬ 
tion  of  (2.1)  is  given  by 


p  =  (7  -  l)p(e  -  i«,u,),  r  =  (e  -  |ui«i)/C,  , 

(2.3) 

which  are  valid  for  perfect  gas,  where  7  is  the  ratio  of 
the  specific  heats  and  is  the  specific  heat  at  con¬ 
stant  volume.  The  components  of  the  viscous  stress 
tensor  (r,,-  are  given  by 


,dui  .  flu,- 


*1 


The  thermal  conductivity  k  and  viscosity  coefficient  p 
are  assumed  to  be  a  function  of  the  temperature  and 
determined  using  Sutherland’s  empirical  relation.  It 
is  assumed  that  A  and  p  are  related  by 


A  = 


2fi 

3 


(2.5) 


find  y  E  T  such  that  'iW  E  W 

(3.1) 

Assuming  Oh  a  classical  triangulation  of  (2  with  the 
nodes  numbered  from  1  to  n  and  Fh  the  boundary  of 
Qh,  we  approximate  the  trial  and  weighting  spaces  T 
and  W  by  their  subspaces  of  finite  dimension  Th  and 
Wfc,  which  respectively,  ate  defined  by 

Th  =  \  Unixj)  I  UH{x,t)  ='^Ui{t)Ni{x) 

I  7=1 

r  "  I 

Wh  .-=  I  IVfc(z)  I  Wh(x)  =  j  (3-2) 


The  left-hand  side  of  equation  (2.1)  constitutes  the  so- 
called  Euler  equations  governing  unsteady  compress¬ 
ible  inviscid  flows. 


where  Nr  is  the  standard  linear  finite  element  shape 
function  associated  with  node  7,  Ui  is  the  value  at 


node  7  and  a/  is  a  constant.  The  Galerkin  finite  ele¬ 
ment  approximation  is  then  given  by 


’  findUh  €  Tfcsuch  that  for  eachiV/(l  <  I  <n) 

'  k  -  fr,  F^mnjNtdTn 

-  /n,  G^m^dClH  -b  &(UH)niNidrH  . 

(3.3) 

The  integrals  appearing  here  are  evaluated  in  the 
standard  finite  dement  form  by  summing  individual 
element  and  boundary  surface  contributions,  the  com¬ 
pact  support  of  the  shape  function  Nj  means  that  the 
equation  can  be  written  as 


F^iUH)njNidrH 


G^mniNrdTH, 


(3.4) 

where  the  summation  extends  over  those  elements  e 
and  boimdary  surfaces  b  that  contain  node  I.  In¬ 
serting  the  assumed  form  for  Uh  in  equation(3.4),  the 
left-hand  side  integral  can  be  evaluated  exactly  to  give 


where  M  denotes  the  finite  element  consistent  mass 
matrix.  For  steady  state  computations,  M  can  be 
replaced  by  the  lumped  (diagonal)  mass  matrix,  de¬ 
noted  by  Ml- 


4.  EDGE-BASED  FINITE  ELEMENT 
SCHEMES  FOR  THE  EULER  EQUATIONS 

It  is  well  known  that  the  discretization  of  the  con¬ 
vective  terms  is  crucial  to  the  successful  numerical  so¬ 
lution  of  the  Navier-Stokes  equations.  Therefore,  we 
start  by  considering  the  solution  of  the  Euler  equa¬ 
tions.  The  extension  to  the  Navier-Stokes  equations 
will  be  discussed  in  the  next  section.  The  integral 
of  the  convective  terms  appearing  on  the  right-hand 
side  of  equation  (3.4)  is  evaluated  approximately  by 
linearly  interpolating  the  flux  in  terms  of  its  nodd 
vdues: 


F’='£^N,F}  ,  F}=F{Ur)  .  (4.1) 

/si 


All  integrals  appearing  in  equation  (3.5)  can  then  be 
evduated  in  closed  form,  leading  to  a  matrix-vector 
product  of  the  form 


Ri  =  DijFi  .  (4.2) 

Here  Djj  are  sparse  matrices,  whose  off-diagond  en¬ 
tries  can  be  identified  as  being  associated  with  the 
edges  of  the  mesh.  Moreover,  it  can  be  shown  that 
for  any  interior  node,  equation  (3.4)  can  be  written 

(4-3) 

SJ 

where  m/  is  the  number  of  edges  connected  to  the 
node  I.  The  coefiScient  Cij  denotes  the  weight  ap¬ 
plied  to  the  average  vdue  of  the  flux  on  the  edge  that 
connects  nodes  I  and  J,  to  obtdn  the  contribution 
made  by  the  edge  to  node  I,  whereas  Cji  denotes  the 
weight  applied  to  the  same  quantity  to  obtain  the  con¬ 
tribution  made  by  the  edge  to  node  J.  These  weights 
are  computed  as 


<^//  =  E 

ee// 


Qe  dNl 
4  dxj 


(4.4) 


where  now  the  summation  extends  only  over  those 
elements  Qe,  which  contain  edge  IJ.  It  can  be  easily 
verified  that  these  weights  possess  the  properties 


mj 

=  0  fordlJ  ,  (4.5) 

IJ 

C{j  =  -Cj,  for  dl  7  and  7  .  (4.6) 

For  notationd  convenience,  we  define  the  vector  Cjj 
by  the  expression 

Clj  =  {C}j,C]j,C]j)  ,  (4.7) 

and  let  Lij  denote  the  modulus  and  Sij  denote  a  unit 
vector  in  the  direction  of  Cu ,  then  equation  (4.3)  can 
be  written  as 


dll 

ij 

where 

F,  =  {F},F},Ff).Su  (4.9) 

Fj={F},F},F^).Su  .  (4.10) 


4.1  Edge-Based  Data  Structure 
The  dternative  procedure  for  obtaining  the  dis¬ 
crete  form  of  the  equations  is  now  apparent.  While 
with  the  element-based  data  structure  information  is 
gathered  from  all  the  nodes  of  each  element,  oper¬ 
ated  on  the  element,  and  then  scattered  back  to  the 
nodes  of  the  element,  the  edge-based  dgorithm  gath¬ 
ers  information  from  dl  the  nodes  of  each  edge,  op¬ 
erates  it  on  the  edge,  and  then  scattes  it  back  to 
the  nodes  of  the  edge.  A  significant  reduction  in 
gather /scatter  costs  and  memory  requirements  can 
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be  realized  by  going  from  an  element-based  to  an 
edge-based  daoa  structure.  For  a  typical  mesh  of 
triangle  elements  with  N  nodes,  the  number  of  ele¬ 
ments,  Nelem  is  about  2N  and  the  number  of  edges, 
Nedge  3N.  An  element-based  data  structure  requires 
2*3*Nelem,  i.e.,  12N  gather/scatter  operations,  while 
its  edge-based  counterpart  needs  2'*‘2*Nedge,  i.e.  12N 
gather/scatter  operations.  Whereas  there  appears  to 
be  little  difference  between  an  edge-based  and  an  el¬ 
ement  based  data  structure  in  2-D,  the  situation  is 
different  in  3-D.  For  a  typical  tetrahedral  mesh  with 
N  nodes,  the  number  of  elements,  Nelem  is  about  5.5N 
and  the  number  of  edges,  Nedge  7N.  An  element- 
based  data  structure  requires  2*4*Nelem,  i.e.,  44N 
gather/scatter  operations,  while  its  edge-based  coun¬ 
terpart  needs  2*2*Nedge,  i.e.,  28N.  Note  that  a  sig¬ 
nificant  gather/scatter  overhead  reduction  is  achieved 
using  an  edg-based  data  structure  in  3D,  thus  leading 
to  a  remarkable  CPU  savings. 

The  memory  overhead  for  an  element-based  data 
structure  calls  for  storing  the  derivatives  of  the  shape 
functions  and  the  volume  of  elements,  and  requires 
a  7*Nelem,  i.e.  14N  vector  in  2-D,  and  a  13*Nelem, 
i.e.  71.5N  vector  in  3-D.  The  memory  overhead  for 
an  edge-based  data  structure  calls  for  storing  the 
first  order  derivatives  and  the  Laplacan  and  needs  a 
3*Nedge,  i.e.  9N  vector  in  2-D,  and  a  4*Nedge,  i.e. 
28N  vector  in  3-D.  Thus  an  edge-based  data  struc¬ 
ture  requires  significantly  less  storage  overhead.  In 
fact,  the  edge-based  data  structure,  defined  by  a  list 
of  edges  with  the  addresses  of  two  nodes  delimiting 
each  edge,  representes  the  minimum  amount  of  infor¬ 
mation  required  to  described  the  unstructured  mesh. 

In  addition,  it  will  be  shown  below  that  this  edge- 
based  data  structure  is  extremely  useful  in  the  process 
of  constructing  different  numerical  schemes,  and  pro¬ 
vides  a  unified  approach  in  which  the  relation  between 
centered  rmd  upwind  schemes  is  clearly  apparent.  It 
is  clear  that  equation  (4.8)  is  nothing  but  a  classic 
Galerkin  finite  element  scheme,  which  is  equivalent 
to  a  central  difference  type  scheme.  By  using  the 
results  of  equation  (4.6),  this  scheme  allows  for  the 
appearance  of  chequerboarding  modes  and  thus  suf¬ 
fers  &om  numerical  instabilities,  unless  some  type  of 
numerical  dissipation  in  the  form  of  artificial  viscos¬ 
ity  is  introduced.  To  construct  stable  schemes  for  the 
Euler  equations,  we  have  to  replace  the  actual  flux 
function  Fu  by  a  consistent  numerical  flux  Tij,  and 
write  the  right-hand-side  in  the  form 


m/ 

=  •  (411) 

IJ 


Then,  by  adopting  different  forms  for  this  numerical 
flux  function,  we  are  able  to  construct  a  number  of  dif¬ 
ferent  algorithms  for  the  Euler  equations  as  described 
below. 


4.2  Godunov^  Scheme  -  Exact  Riemann  Solver 
A  stable  scheme  can  be  constructed  by  using  the 
exact  Riemann  solver  [10].  This  would  imply  replac¬ 
ing  (4.11)  by 

mi 

RHS{I)  =  '£LaF{U,^^  ^4.12) 

IJ 

where  Ufj  denotes  the  local  exact  soli  ..>n  of  u^e  Rie¬ 
mann  problem  to  the  Euler  equation  and  can  be  ex¬ 
pressed  as 

U^  =  RieiU,,Ur)  (4.13) 

where  we  have  set 

Ur  =  Ui,  U,  =  Uj  .  (4.14) 

This  is  the  first  order  Godunov  scheme.  A  scheme 
of  higher  order  accuracy  can  be  achieved  by  a  better 
approximation  to  Ur  and  Ui,  e.g.,  via  reconstruction 
process  and  monotone  limiting.  The  major  disadvan¬ 
tage  of  Godunov’s  approach  is  the  extensive  computa¬ 
tional  work  introduced  through  the  Riemann  solver. 

4.3  Roe  Scheme  -  Approximate  Riemann  Solver 
A  first  simplification  can  be  achieved  by  replac¬ 
ing  the  computationally  costly  exact  Riemann  solver 
by  an  approximate  Riemann  solver.  A  variety  of  pos¬ 
sibilities  can  be  defined;  here  we  consider  one  of  the 
most  popular  approximate  Riemann  solvers,  namely 
the  flux-difference  splitting  of  Roe  [11]: 


=  Fi  +  Fj-  1  Aij  I  (Uj  -  Ur)  .  (4.15) 

Here  |  Atj  j  denotes  the  standard  Roe  matrix  eval¬ 
uated  in  the  direction  Srj-  It  can  be  shown  that 
this  scheme  is  equivalent  to  the  first  order  finite  vol¬ 
ume  upwind  cell-vertex  scheme  based  on  a  dual  mesh. 
Many  different  ways  exist  to  achieve  higher  order  ac¬ 
curacy.  In  the  present  study,  a  scheme  of  higher  or¬ 
der  accuracy  is  obtained  by  using  upwind-biased  in¬ 
terpolations  of  the  solution  U  via  the  MUSCL  ap¬ 
proach  [12].  This  leads  to  the  flux  function 

Fu  =  F+  +  FJ-  I  A(U+,  UJ)  1  {UJ  -  Ut)  (4.16) 

where 

n  =  mt),  FJ  =  F{UJ)  .  (4.17) 

The  upwind-biased  interpolations  for  Uf  and  UJ  are 
defined  by 

i/+  =  U,  +  i[(l  -  *)A7  (H-  k){Uj  -  Ur)]  (4.18) 

UJ  =  Uj-  i[(l  -  fc)Aj  -1-  (1-1-  k){Uj  -  Ui)]  (4.19) 
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'  where  the  forward  and  backward  difference  operators 
are  given  by 

a;  =  Ui-  Ui.i  =  2(vtr)/  •  i"  -  {Uj  -  u,)  (4.20) 

AJ  =  f/j+i  -Uj  =  2{VU)j  ■  1"  -  {Uj  -  Ui)  (4.21) 

where  =x/  —  x/  is  the  length  vector  of  this  edge. 

The  parameter  i,  which  can  be  chosen  to  control 
the  degree  of  approximation,  was  set  to  k  =  1/3  in 
all  calculations  presented  here.  This  correponds  to  a 
third-order  upwind-biased  scheme  [12].  With  higher 
order  spatial  accuracy,  spurious  oscillations  in  the 
vicinity  of  shock  waves  are  expected  to  occur.  Some 
form  of  limiting  is  usually  required  to  eliminate  these 
numerical  oscillations  of  the  solution  and  to  provide 
some  kind  of  monotonicity  property.  The  flux  lim¬ 
iter  modifies  the  upwind-biased  int^olation  f//  and 
C/j  and  the  equations  (4.18)  and  (4.19)  are  replaced, 
respectively  by 

Ut  =  Ut  +  ^[(1  -  fcs,)A7  +  (1  -I-  k8i)iUj  -  Ui)] 

(4.22) 

uy  =  Uj-  ^[(1  -  k3j)A^  +  (1  +  ksj)(Uj  -  Ui)] 

(4.23) 

where  s  is  the  flux  limiter.  Both  the  Van  Albada  lim¬ 
iter  and  Miiunod  limiter  were  employed  in  this  study. 
Three  options  exist  concerning  the  choice  of  interpola¬ 
tion  variables:  conservative  variables,  primitive  vari¬ 
ables,  and  characteristic  variables.  Using  limiters  on 
characteristic  variables  seems  to  give  the  best  results, 
but  sacrifices  some  computational  ef&ciency. 

4.4.  Scalar  Limited  Dissipation 
A  further  possible  simplification  can  be  made  by 
replacing  the  Roe  matrix  by  its  spectral  radius.  This 
le^  to  the  choice 


=  Fi  +  Fj-  I  Xjj  I  iUj  -  Ui)  (4.24) 
for  the  numerical  flux  function,  where 

I  Xu  1=1  ujj  •  Sjj  I  +cij  (^’SS) 

and  “/J  and  cij  denote  edge  values,  computed  as 
nodal  average,  of  the  fluid  velocity  and  speed  of 
sound  respectively.  This  can  be  considered  as  a  cen¬ 
tered  difference  scheme  plus  a  second  order  dissi¬ 
pation  operator,  leading  to  a  first  order,  monotone 
scheme.  A  higher-order  scheme  would  be  obtained  by 
a  better  approximation  to  the  ‘right’  and  ‘left’  states 
of  the  ‘Riemann  problem’,  which  have  been  set  to 
Ur  =  Ui,Ui  =  Uj.  This  would  reduce  the  difference 
between  Ui,Uj,  decreasing  in  turn  the  dissipation. 
As  before,  limiting  would  provide  automatic  cut-off 
values  for  the  possible  range  of  Uj,Uj  that  lead  to 
monotonic  solutions. 

4.5.  Scalar  Dissipation  With  Pressure  Sensors 
Another  option  to  reduce  the  magnitude  of  — 

1  is  to  apply  a  blended  second-  and  fourth-order 


damping  for  the  dissipation  [Ij.  This  is  borne  by  the 
observation  that  even  for  smooth  problems,  central 
difference  schemes  still  require  fourth  order  damping 
for  stability.  A  scheme  of  this  type  may  be  written 
as: 

Tij  =  Fi  +  Fj 


-\hj\ 


Ui-Uj-¥  |l"(VC^/  -I-  VUj) 


(4.26) 


where  0  <  /?  <  1  denotes  a  pressure  sensor  function 
of  the  form  [1] 


PI -Pi-  +  ^PJ)  u  vi\ 

^  |p/-p/H-|0.51"(Vp/-t-Vpj)l  • 

We  have  experimented  with  several  combinations  of 
this  sensor  function.  The  one  we  favor  at  the  present 
time  is  computed  in  two  passes  over  the  mesh.  In  the 
first  pass,  the  highest  of  the  edge-based  ^-values  given 
by  eqnuation  (4.27)  is  kept  for  each  point.  In  a  second 
pass,  the  highest  v^ue  of  the  two  points  belonging  to 
an  edge  is  kept  as  the  final  /?- value.  For  0  =  0, 1, 
second-  and  fourth-  order  damping  operators  are  ob¬ 
tained  respectively.  We  remark  that  although  this 
discretisation  of  the  Euler  fluxes  looks  like  a  blend 
of  second-  and  fourth-order  dissipation,  it  has  no  ad¬ 
justable  parameters. 

4.6.  Scalar  Dissipation  Without  Gradients 
The  scalar  dissipation  operator  presented  above 
still  requites  the  ev^uation  of  gradiaits.  This  can 
be  quite  costly  for  Euler  simulations.  An  alterna¬ 
tive  is  the  simplify  the  combination  of  second-  and 
fourth-order  damping  operators  by  writing  out  explic¬ 
itly  these  operators: 

d2^Xu{l-0)[Ui-Uj]  , 


dA  =  Xu0  [c^j  -Uj-^  • 

Performing  a  Thylor-series  expansion  in  the  direction 
of  the  edge,  we  have 


1"  P  \SFU  8^U  1 

Ui-Uj^\{VUi+VUj)  «  4  [9^  I.  ~  l/J  . 

This  suggests  the  following  simplification,  which  ne¬ 
glects  the  off-diagonal  terms  of  the  tensor  of  second 
derivatives: 


f*  \d^U  ,  d^U  ,]  r„2,,  _2„  1 


and  leads  to  the  familiar  blend  second-  and  fourth- 
order  damping  operators 
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=  F/  +  Fj-  I  A//  1(1-/?)  [Ui  -  Uj] 

-A/j^^  [V^Uj  -  V^t/j]  (4.28) 

4.7  Taylor- Galerldn  scheme 

Due  to  their  importance  for  transient  calcula¬ 
tions,  it  is  worthwhile  to  consider  possible  edge-based 
Taylor-Galerkin  schemes.  The  essential  feature  of  any 
Taylor-Galerkin  scheme  is  the  combination  of  time 
amd  space  discretizations,  leading  to  second-order  ac¬ 
curacy  in  both  time  and  space.  An  edge-based  two 
step  Taylor-Galerkin  scheme  can  be  readily  obtained 


by  adopting  the  numerical  flux 

(4.29) 

where 

(4.30) 

)lj  is  computed  on  each  edge  and  given  by 


The  major  advantage  of  this  scheme  is  its  speed,  since 
there  is  no  requirement  of  gradient  computations  and 
limiting  procedures.  An  explicit  numerical  dissipation 
in  the  form  of  Lapidus  viscosity  is  needed  to  model 
flows  involving  di^ntinuties.  The  Taylor-Galerkin 
scheme  alone  is  of  little  use  practically.  However  it 
provides  a  useful  base  scheme  for  the  flux-corrected 
transport  scheme  presented  below. 

4.8  Flux-Corrected  Transport  Scheme 

The  idea  behind  FCT  is  to  combine  a  high-order 
scheme  with  a  low-order  scheme  in  such  a  way  that 
the  high-order  scheme  is  employed  in  smooth  regions 
of  the  flow,  whereas  the  low-order  scheme  is  used  near 
discontinuities  in  a  conservative  way,  in  an  attempt  to 
yield  a  monotonic  solution.  The  implementation  of 
an  edge-based  FCT  scheme  is  exactly  the  same  as  its 
element-based  counterpart  [13].  However,  the  use  of 
an  edge-based  data  structure  makes  the  implemen¬ 
tation  more  efficient,  which  is  especially  attractive 
for  three  dimensional  problems.  As  the  high-order 
scheme,  we  employ  the  edge-based  two  step  Taylor- 
Galerkin  scheme  with  consistent  mass  matrix.  The 
converged  solution  can  be  recast  into  the  following 
form: 


MLAUh  =  R  +  (Ml- Mc)AU''  .  (4.32) 


i.e.,  lumped  mass-matrix  plus  mass  diffusion.  Sub¬ 
tracting  (4.33)  from  (4.32)  yields  the  antidiifusive 
edge  contributions 

(At/'*  -  At/')  =  MI\Ml  -  Mc){  .7"  -  vt/'*)  . 

(4.34) 

This  avoids  any  need  for  physical  flux  ^coi  stations 

and  leads  to  a  very  fast  overall  s'  ne.  though 

FEM-FCT  is  often  criticized  as  n  aav:  \  strict 

mathematical  background,  our  ex  lenc-  xs  been 

that  it  gives  excellent  resolution  f  both  jcks  and 
contact  discontinuities  for  the  si  olation  of  strongly 
unsteady  compressible  flows. 

5.  EXTENSION  TO 
THE  NAVIEB^STOKES  EQUATIONS 

In  this  section,  we  extend  the  solution  algorithms 
described  above  for  the  Euler  equations  to  the  solu¬ 
tion  of  the  Navier-Stokes  equations.  The  numerical 
integration  of  the  viscous  terms  is  carried  out  in  a 
centered  way.  Note  that  the  viscous  fluxes  are  func¬ 
tions  of  unknowns  and  their  first  derivatives,  which 
can  be  expressed  as 


(5.1) 


The  following  are  three  possible  ways  to  evaluate  the 
integrals  that  involve  these  viscous  fluxes. 

(a)  Standard  Finite  Element  Approach 

In  this  approach,  the  integral  involving  the  vis¬ 
cous  fluxes  is  equated  numerically  after  inserting  the 
assumed  form  for  U  into  the  viscous  flux  functions. 
Note  that  over  each  element,  dUfdx^  will  be  constant 
since  U  is  assumed  to  have  a  linear  variation.  Thus 
for  an  element  e  with  nodes  I,  J,  K  and  L,  the  inte^al 
is  obtained,  using  a  one  point  numerical  integration 
rule,  as 

U  CPi\{Ui+Uj+UK+Ui.),^  u) .  (5.2) 


(b)  Mixed  Formulation 

An  alternative  approach  is  to  evaluate  the  vis¬ 
cous  flux  contributions  by  making  using  of  the  nodal 
gradients  of  the  unknown  vector  U.  In  this  case,  nodal 
values  of  the  viscous  fluxes  can  be  directly  evaluated 


as 


G>  =  G^(f/.£). 


Now,  the  viscous  fluxes  can  be  interpolated  linearly 
over  each  element  and  considered  jointly  with  the  in- 
viscid  fluxes.  In  this  case,  of  course,  the  edge  data 
structure  can  be  readily  employed  to  give  the  approx¬ 
imation 


The  low  order  scheme  used  is  simply 

Ml^Ui  =  R  +  C4{Mc-Ml)U^  , 


(4.33) 
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>  It  is  noted  here,  that  the  discretization  of  the  viscous 
terms  using  the  standard  finite  element  approach  only 
uses  information  from  those  points  which  are  directly 
connected  to  the  points  being  considered.  On  the 
other  hand,  the  mixed  formulation  involves  informa¬ 
tion  from  two  layers  of  points  surrounding  the  point 
under  consideration.  Despite  their  inconsistency,  the 
numerical  experience  indicates  that  the  results  from 
two  approaches  are  virtually  identical. 

(c)  Edge-based  Finite  Element  Approach. 

Note  that  the  integrals  of  viscous  flux  terms  in¬ 
volves  the  evaluation  of  the  following  terms 


f  .  dg  dNi 
'Hk  ^  dxk  dzj 


dQn  . 


It  is  easily  shown  that  this  integral  can  be  written  as 


dg  dNi 
dxk  dxj 


mi 

dQH  =  Yi 


ij 


where  Du  are  computed  as 


(5.5) 


teiJ 


Cle  dNi 
4  dxt 


dNi 

dxj 


(5.6) 


This  approach  provides  a  way  to  evaduate  the  inte¬ 
gral  of  the  viscous  fluxes  that  is  consistent  with  the 
standard  finite  element  approach.  The  computational 
effort  for  this  approach  is  very  small.  The  major  dis¬ 
advantage  of  this  approach  is  that  overhead  storage 
requirements  become  significant. 


6.  TEMPORAL  DISCRETIZATION 

Equation(4.8)  represents  the  time  evalution  of 
the  unknown  vector  Ui{t)  at  node  I  of  the  grid.  As¬ 
suming  that  the  nodal  v^ues  Uj  are  known  at  time 
i„,  the  solution  is  advanced  over  a  time  step  At, 
to  time  by  an  explicit  multi-stage  Runge-Kutta 
time-stepping  scheme  given  by 

u\°'>  =  uy 


[/(P)  =  f;W-apAf(A/i)7'R/(£/('’-‘))  P= 


right-hsmd  side  might  be  evaluted  only  for  the  first 
stage. 


7.  ADAPTIVE  REFINEMENT 

A  very  attractive  feature  of  unstructured  grids 
is  the  ease  with  which  they  incorporate  adaptive  re¬ 
finement.  The  addition  of  further  degrees  of  free¬ 
dom  does  not  destroy  any  previous  structure.  Thus, 
the  flow  solver  requires  no  further  modification  when 
operating  on  an  adapted  grid.  For  many  practical 
problems,  the  regions  that  need  to  be  refined  are  ex¬ 
tremely  small  as  compared  to  the  overall  domain.  On 
the  other  hand,  the  spatial  location  of  these  regions 
where  small  elements  are  required  are  typically  un¬ 
known,  and  may  vary  in  time.  It  is  thus  not  surpris¬ 
ing  that  the  use  of  Captive  refinement  typically  ac¬ 
crues  savings  in  storage  and  CPU  requirements,  which 
range  between  10-100  as  compared  to  an  overall  fine 
mesh.  Any  adaptive  refinement  scheme  consists  of 
three  different  stages:  determining  what  we  want  to 
achieve  by  refining  the  grid,  developing  an  error  in¬ 
dicator/estimator  to  detect  the  regions  to  be  refined, 
and  a  refinement  strategy,  such  as  movement,  enrich¬ 
ment  or  remeshing.  In  this  study,  an  h-refinement 
scheme  has  been  incorporated  into  the  flow  solvers  in 
order  to  enhance  the  solution  accuracy  and  efficiency. 
Further  detail  and  description  about  this  adaptive 
scheme  can  be  found  in  [14-15]. 

8.  NUMERICAL  EXAMPLES 

The  results  obtained  by  the  centered  difference 
schemes  are  not  included  in  this  paper,  since  such 
results  can  be  found  in  [1]  for  ste^y  state  compu¬ 
tations,  and  more  recently  in  [16]  for  the  simulation 
of  high  speed  trains  through  tunnels  using  an  ALE 
formulation  and  adaptive  remeshing  techniques.  As 
the  solutions  obtain^  by  the  edge-based  FEM-FCT 
scheme  are  almost  identical  to  those  obtained  by  its 
element-based  counterpart,  only  a  few  test  cases  are 
presented  for  the  edge-based  FEM-FCT  scheme.  All 
the  results  to  be  presented  are  obtained  by  edge-based 
upwind  finite  element  scheme,  where  the  Van  Albada 
limiter  based  on  primitive  variables  is  used  for  the 
steady  state  computations  and  the  Minmod  limiter 
based  on  characteristic  variables  is  used  for  transient 
flow  calculations.  For  the  purpose  of  comparison  , 
the  simulations  for  test  cases  3-4  were  also  perfomed 
using  the  edge-based  FEM-FCT  scheme  to  show  its 
excellent  resolution  for  both  shock  waves  and  conttict 
discontinuities  for  the  simulation  of  strongly  unsteady 
compressible  flows. 


with  the  parameters  Op  assigned  appropriate  values. 
The  scheme  is  second  order  accurate  in  time.  For 
steady  state  computations,  implicit  residual  smooth¬ 
ing  and  local  time-stepping  are  used  to  accelerate 
convergence  to  steady  state.  Residual  smoothing  2d- 
lows  the  use  of  larger  CFL  numbers  than  the  one  dic¬ 
tated  by  the  stability  of  the  original  scheme.  For  the 
centered-difference  scheme,  in  the  interests  of  compu¬ 
tational  efficiency,  the  diffusion  contribution  to  the 


Test  Case  1.  NACA0012  Airfoil 

The  problem  under  consideration  is  transonic 
flow  around  a  NACA0012  airfoil  with  a  freestream 
Mach  number  of  0.85  and  an  angle  of  attack  of  1  de¬ 
gree  -  a  classical  test  problem  for  Euler  solvers.  The 
grid  adaption  scheme  was  used.  The  initial  mesh  con¬ 
taining  1,311  elements  and  719  points  and  the  final 
refined  mesh  consisting  of  6,397  elements  and  3,274 
points  after  three  levels  of  refinement  are  shown  in 
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Figures  la  and  lb,  respectively.  The  computed  pres¬ 
sure  contours  in  the  flow  fields  on  the  initial  and  final 
adapted  mesh  are  displayed  in  Figures  Ic  and  Id,  re¬ 
spectively.  Figures  le  and  If  show,  respectively,  the 
comparison  of  pressure  coefficients  and  entropy  on  the 
airfoil  between  initial  mesh  and  final  refined  mesh. 
The  refinement  of  regions  with  strong  gradients  such 
la  shocks,  leading  edge,  and  trailing  edge  is  well  pre¬ 
sented  and  the  considerable  improvement  in  the  solu¬ 
tion  for  these  regions  after  refinement  is  clearly  appar¬ 
ent.  The  production  of  numerical  entropy  in  the  vicin¬ 
ity  of  the  stagnation  point  is  dramatically  decreased 
after  refinement.  The  pressure  coefficient  distribution 
on  both  initial  and  final  adapted  mesh  indicates  that 
there  is  only  one  grid  point  within  the  shock  structure 
and  demonstrates  the  sharp  shock  capturing  ability  of 
Roe’s  approximate  Riemann  solver  for  the  solution  of 
steady  problems. 

Test  Case  2.  Half  Cylinder 

In  this  example,  we  consider  hypersonic  flow 
around  a  half  cylinder  with  a  freestream  Mach  num¬ 
ber  of  10.  The  particular  difficulty  is  due  to  the  large 
Mach  number  and  a  quasi-rarefaction  zone  behind  the 
cylinder.  The  final  refined  mesh  after  three  levels  of 
refinement,  shown  in  Figure  2a,  contains  20,178  ele¬ 
ments  and  10,277  points.  The  computed  Mach  num¬ 
ber  contours  in  the  flow  fields  are  depicted  in  Figure 
2b.  The  pressure  coefficient  distribution  and  entropy 
on  the  surface  are  shown  in  Figure  2c  and  2d,  respec¬ 
tively.  For  this  computation,  the  choice  of  Mach  num¬ 
ber  as  error  indicator  provides  the  best  refinement  in 
the  wake  of  blunt  bodies. 

Test  Case  3.  Shock  Tube  Problem  or  Riemann  Problem 

The  shock  tube  problem  constitutes  a  partic¬ 
ularly  interesting  and  difficult  test  case,  since  it 
presents  an  exact  solution  to  the  full  system  of  one¬ 
dimensional  Euler  equations  containing  simultane¬ 
ously  a  shock  wave,  a  contact  discontinuity,  and  an 
expansion  fan.  The  initial  conditions  in  the  present 
computation  are  the  following: 

p  =  1.000,ut  =  0,p  =  1.0,  0.0<x<50.0 

p  =  0.125, u/j  =  0,p  =  0.1,  50.  <  z  <  100. 

Figure  3  shows  the  computational  mesh  and  the 
results  for  the  edge-based  FEM-FCT  scheme  and 
edge-based  upwind  finite  element  scheme.  This  is 
a  2D  simulation  of  a  ID  problem.  The  mesh  con¬ 
sists  of  101  points  in  the  X-direction  and  3  points  in 
the  Y-direction.  Both  schemes  produced  very  similar 
results,  with  excellent  resolution  for  both  shock  and 
contact  discontinuity. 

Test  Case  4.  Shock  Diffraction  in  a  Baffled  T\ibe 

The  problem  under  consideration  is  shown  in  Fig¬ 
ure  4.  A  weak  shock  {M,  =  1.31)  propagates  inside 
a  baffled  tube.  We  selected  this  problem  to  show 
that  the  FCT  scheme  together  with  the  classic  h- 
enrichement /coarsening  grid  adaption  scheme  could 
produce  excellent  results  for  the  efficient  simulation 
of  strongly  unsteady  flows,  without  deteriorating  the 
accuracy  of  the  solutions.  The  number  of  refinement 


levels  allowed  in  this  case  is  5  zuid  the  gnd  is  mod¬ 
ified  every  7  timesteps.  Density  was  chosen  as  the 
key  variable  for  the  error  indicator.  Figure  4a  shows 
the  adapted  mesh  at  60  ps  and  the  experimental  solu¬ 
tion  is  depicted  in  Figure  4b.  Figures  4c  and  4d  show 
the  computed  density  contours  obtained  by  the  edge- 
based  FEM-FCT  and  the  edge-based  upw^  d  finite 
element  ^  rheme,  respectively.  Comparisor  f  these 
results  demonstrates  that  the  edge-based  ’  .M-FCT 
producea  better  solutions  than  the  edge-  ^sed  up¬ 
wind  finite  element  scheme,  while  runni  approxi¬ 
mately  3  times  faster. 

Test  Case  5.  ONERA  M6  Wing 

In  this  test  case,  we  consider  a  transonic  flow  over 
the  ONERA  M6  wing  geometry.  The  M6  wing  has  a 
leading  edge  sweep  angle  of  30  degree,  an  aspect  of 
3.8,  and  a  taper  ratio  of  0.562.  The  airfoil  section 
of  the  wing  is  the  ONERA  "D”  airfoil,  which  is  a 
10%  maxinium  thickness-to-chord  ratio  conventional 
section.  The  flow  solutions  are  presented  at  a  Mach 
number  of  0.84  and  an  angle  of  attack  of  3.06.  The 
mesh,  which  contains  179,106  grid  points,  951,179  el¬ 
ements,  and  34,013  boundary  points,  is  depicted  in 
Figure  5a.  Figures  5b  and  5c  show  pressure  contours 
on  the  upper  wing  surface  and  lower  surface,  respec¬ 
tively.  The  upper  surface  contours  clearly  show  the 
sharply  captured  lambda-type  shock  structure  formed 
by  the  two  inboard  shock  waves,  which  merge  together 
near  80%  semispan  to  form  the  single  strong  shock 
wave  in  the  outboard  region  of  the  wing.  The  com¬ 
puted  pressure  coefficient  distributions  are  compared 
with  experimented  data  [17]  at  six  spanwise  stations 
in  Figure  5d.  The  results  obtained  compare  closely 
with  experimental  data,  except  at  the  root  stations, 
due  to  lack  of  viscous  effects. 

Test  Case  6.  Boeing  747  Aircraft 

The  sixth  test  case  is  performed  on  a  com¬ 
plete  Boeing  747  aircraft.  The  747  configuration  in¬ 
cludes  the  fuselage,  the  wing,  horizontal  and  vertical 
tails,  underwing  pylons,  and  flow-through  engine  na¬ 
celles.  The  mesh,  which  contains  162,440  grid  points, 
901,275  elements  and  18,454  boundary  points  for  the 
half-span  airplane,  is  shown  in  Figure  6a.  Solutions 
were  computed  for  the  aircraft  at  a  free  stream  of 
Mach  number  of  0.84  and  an  angle  of  attack  of  2.73 
degrees.  The  computations  were  performed  using  the 
three-stage  Runge-Kutta  time-stepping  scheme  with 
local  time  stepping  and  implicit  residual  smoothing. 
The  solution  was  advanced  100  time-steps  using  a 
CFL  number  of  0.35  and  1000  time-steps  with  a  CFL 
number  of  4.0,  to  converge  the  solution  to  engineering 
accuracy(  a  decrease  of  a  three  order-of-magnitude  in 
the  L2  norm  of  the  density  residual).  The  computed 
pressure  contours  on  the  surface  of  the  airplane  are 
shown  in  Figure  6b. 

Test  Case  7.  Viscous  Flow  Past  a  Flat  Plate 

This  test  case  involves  a  laminar  flow  past  a  flat 
plate  at  a  Mach  number  of  0.5  and  a  chord  Reynolds 
number  of  10,000.  This  computation  was  performed 
to  validate  the  Navier-Stokes  code  by  comparing  the 
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■  results  with  the  exact  Blasius  solutions.  The  mesh 
used  in  the  computation  is  shown  in  Figure  7a.  It  con¬ 
tains  2,604  elements,  1,376  points,  and  146  boundary 
points.  The  computed  Mach  nun^er  contours  in  the 
flow  field  are  depicted  in  figure  7b,  where  the  devel¬ 
opment  of  a  boundary  layer  can  be  clearly  observed. 
Figure  7c  shows  the  comparison  of  the  Blasius  veloc¬ 
ity  profile  and  the  computed  velocity  profiles  as  scaled 
by  the  Blasius  similarity  law  at  different  chord  length 
downstream  of  the  lea^g  edge.  The  computed  re¬ 
sults  indicate  that  the  similarity  solution  for  a  flat 
plate  boundary  layer  is  correctly  obtained  and  the  so¬ 
lution  agrees  well  with  the  Blasius  solution.  Finally, 
Figure  7d  shows  the  comparison  of  velocity  profiles 
obtained  using  two  different  discretization  approachs 
for  viscous  terms.  It  is  observed  that  these  two  types 
of  discretization  give  practically  identical  results. 

7.  CONCLUSIONS 

We  have  given  an  overview  of  the  recent  devel¬ 
opment  of  some  high  order  accuracy  finite  element 
schemes  to  the  solutions  of  the  compressible  Euler  and 
Navier-Stokes  equations  on  unstructured  grids.  These 
schemes  are  bas^  on  an  edge-based  data  structure,  as 
opposed  to  the  more  traditionel  element-based  data 
structure.  The  advantages  of  such  data  structure 
might  be  summarized  as  a)  better  performance  of 
the  numerical  schemes  in  terms  of  computational  ef¬ 
ficiency  and  memory  requirements,  b)  easy  construc¬ 
tion  of  different  numerical  schemes,  ranging  from  the 
Godunov  scheme  to  the  centered  scheme  with  blended 
dissipation,  and  c)  easy  implementation  of  numerical 
schemes  for  both  2D  and  3D  problems,  due  to  1-D 
like  data  structure. 
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Fig.la  Initial  Mesh  Used  for  Computing  Flow 
past  NACA0012  Airfoil,  nelem=1311,  npoin=719 


Fig.lb  Final  Adapted  Mesh  Used  for  Computing  Flow 
past  NACA00I2  Airfoil,  nelem=6397,  npoin=3274 


Fig.lc  Computed  Pressure  Contours 
on  Initial  Mesh,  Mq,  =  0.85,  o  =  1 


Fig.ld  Computed  Pressure  Contours  on 
the  Final  Adapted  Mesh,  Moo  =  0.85,  a  = 


Fig.le  Comparison  of  C^  Distribution  Obtained 
Using  the  Initial  and  Final  Meshes 


Fig.  If  Comparison  of  Entropy  Distribution 
Obtained  Using  the  Initial  and  Final  Meshes 


Fig.4a  Adapted  Mesh  at  60  ^ls  for  Computing 
Shock  diffraction  in  a  Baffled  liibe 
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Fig.7a  Mesh  Used  for  Computing  Viscous  Flow 
past  a  Flat  Plate;  nelem=2,604,  npoinsl,376 


Fig.7b  Computed  Mach  Number  Contours  past 
a  Flat  Plate;  M«os0.5,  Res:10,000 
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Fig.7c  Streamwise  Velocity  Profiles  past 
a  Flat  Plate;  Meo=:0.5,  Ile=10,000 


Fig.7d  Comparison  of  Velocity  Profiles  Obtained 
Using  Different  Approachs  to  Model  Viscous  Terms 


17 


AM  A  A 


AIAA  93-2933 

Numerical  Solution  of  the  Euler 
Equations  for  Complex  Aerodynamic 
Configurations  Using  an  Edge-Based 
Finite  Element  Scheme 

Hong  Luo*, Joseph  D.  Baum*,  and  Rainald 
L6hner** 

*  Science  Applications  International  Corporation, 
McLean,  VA 

**  George  Washington  University, 

Washington,  D.C. 


AIAA  24th 

Fluid  Dynamics  Conference 

July  6-9, 1993  /  Orlando,  FL 


For  permission  to  copy  or  republish,  contact  the  American  institute  of  Aeronautics  and  Astronautics 
370  L'Enfant  Promenade,  S.W.,  Washington,  D.C.  20024 


NUMERICAL  SOLUTION  OF  THE  EULER  EQUATIONS 
FOR  COMPLEX  AERODYNAMIC  CONFIGURATIONS 
USING  AN  EDGE-BASED  FINITE  ELEMENT  SCHEME 


Hong  Luo  and  Joseph  D.  Baum 
Science  Applications  International  Corporation 
1710  Goodridge  Drive,  MS  2-3-1 
McLean,  VA  22102,  USA 

and 

Rainald  Lohner 

CMEE,  School  of  Engineering  and  Applied  Science 
The  George  Washington  University,  Washington,  D.C.  20052,  USA 


ABSTRACT 

This  paper  describes  the  vdevelopment,  validation 
and  application  of  a  new  finite  element  scheme  for 
the  solution  of  the  compressible  Euler  equations  on 
unstructured  grids.  The  implementation  of  the  nu¬ 
merical  scheme  is  based  on  an  edge-based  data  struc¬ 
ture,  as  opposed  to  a  more  traditional  element-based 
data  structure.  The  use  of  this  edge-based  data  struc¬ 
ture  not  only  improves  the  efficiency  of  the  algorithm, 
but  also  enables  a  straightforward  implementation  of 
upwind  schemes  in  the  context  of  finite  element  meth¬ 
ods.  The  algorithm  has  been  tested  and  validated  on 
some  well  documented  configurations.  A  flow  solution 
about  a  complete  F-18  fighter  is  shown  to  demon¬ 
strate  the  accuracy  and  robustness  of  the  proposed 
algorithm. 

1.  INTRODUCTION 

In  recent  years,  significant  progress  has  been 
made  in  the  development  of  numerical  algorithms  for 
the  solution  of  the  compressible  Euler  and  Navier- 
Stokes  equations.  The  use  of  unstructured  meshes  for 
computational  fluid  dynamics  problems  has  become 
widespread  due  to  their  ability  to  discretize  arbitrar¬ 
ily  complex  geometries  and  due  to  the  ease  of  adap¬ 
tion  in  enhancing  the  solution  accuracy  and  efficiency 
through  the  use  of  adaptive  refinement  techniques. 
However,  any  unstructured  algorithm  requires  the 
storage  of  the  mesh  connectivity,  which  implies  the 
increase  of  computer  memory  and  the  use  of  indi¬ 
rect  addressing  to  retrieve  nearest  neighbor  informa¬ 
tion.  These  requirements,  in  turn,  mean  that  any 
numerical  algorithm  will  run  slower  on  an  unstruc¬ 
tured  grid  than  on  a  structured  grid.  In  order  to 
reduce  indirect  addressing,  new  edge-based  finite  ele¬ 
ment  8chemes([l]-[4])  have  been  recently  introduced. 
In  addition,  even  more  sophisticated  data  structures 
such  as  stars,  super  edges,  and  chains  were  recently 
developed  by  Lohner[5].  The  use  of  edge-based  data 
structure  has  shown  to  result  in  remarkable  compu¬ 
tational  savings  for  three  dimensional  problems. 

In  the  last  few  years,  extensive  research  has  been 


done  on  upwind  type  algorithms  for  the  solution  of 
the  Euler  and  Navier-Stokes  equations  on  unstruc¬ 
tured  meshes([6]-[9]).  A  significant  advantage  of  up¬ 
wind  discretization  u  that  it  is  naturally  dissipative, 
in  contrast  with  central-difference  discretizations,  and 
consequently  does  not  require  any  problem-dependent 
parameters  to  adjust.  So  far,  all  upwind  schemes 
implemented  as  either  node-centered  or  cell-centered 
discretizations  on  unstructured  meshes  have  used  the 
finite  volume  approach  where  the  control  volume  must 
be  constructed  first.  In  terms  of  computational  effi¬ 
ciency,  node-centered  schemes  are  preferable  to  their 
cell-center  counterparts.  In  the  node-centered  ap- 
ptoach([5],[&]),  the  control  volume  is  typically  taken 
to  be  part  of  the  neighboring  cells  that  have  a  ver¬ 
tex  at  that  node.  In  two  dimensions,  the  part  of  the 
cells  taken  is  determined  by  connecting  the  centroid 
of  the  cell  and  the  midpoints  of  the  two  edges  that 
share  the  node.  In  3-0,  the  part  of  the  cells  taken 
is  determined  by  a  surface  constructed  in  a  similar 
way.  However,  this  is  somewhat  complicated  geomet¬ 
rically  to  do  in  three  dimensions.  The  switching  from 
element  to  edge-based  data  structure  renders  the  im¬ 
plementation  of  upwind  schemes  trivial  and  straight¬ 
forward  in  the  context  of  the  finite  element  approach; 
this  is  especially  attractive  for  three  dimensional  ap¬ 
plication,  since  there  is  no  need  to  construct  control 
volumes  explicitly  and  geometrically. 

The  authors  have  recently  developed  some  high 
sccurru:y  schemes  for  the  solution  of  the  Euler  and 
Navier-Stokes  equations  on  unstructured  grids  by  us¬ 
ing  an  edge-based  data  structure^!].  This  paper  de¬ 
scribes  the  development,  validation,  and  application 
of  an  upwind  finite  element  algorithm  to  the  simula¬ 
tion  of  three  dimensional  compressible  flows  around 
complex  aerodynamic  configurations.  In  this  scheme, 
the  spatial  discretization  is  accomplished  by  an  edge- 
based  finite  element  formulation  using  Roe’s  flux- 
difference  splitting.  A  MUSCL  approi^  is  used  to 
achieve  higher-order  accuracy.  A  characteristic  anal¬ 
ysis  based  on  the  introduction  of  Riemann  invariants 
for  one-dimensional  flow  normal  to  the  boundary  is 
used  to  treat  boundary  conditions.  Solutions  are  ad- 
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vanced  in  time  by  a  multi-stage  Runge-Kutta  time¬ 
stepping  scheme.  Convergence  is  accelerated  using 
local  time-stepping  and  implicit  residual  smoothing. 
The  algorithm  has  been  tested  and  validated  on  some 
well  documented  configurations.  A  solution  of  the 
flow  around  a  complete  F-18  fighter  is  presented  to 
demonstrate  the  accurzury  and  robustness  of  the  pro¬ 
posed  algorithm. 

2.  GOVERNING  EQUATIONS 

The  Euler  equations  governing  unsteady  com¬ 
pressible  inviscid  flows  can  be  expressed  in  the  con¬ 
servative  form  as 


dt  dxj 


=  0, 


(2.1) 


where  the  summation  convention  has  been  employed 
and 


(  P  \  (  P^i  \ 

U  =  \  pui  I  ,F^  =  \  puiUj  +p6ij  I  (2.2) 

\  pe  /  \  Ui(pe  +  p)  J 


where  Nj  is  the  standard  linear  finite  element  shape 
function  associated  with  node  I,  Uj  is  the  value  at 
node  /,  find  a/  is  a  constant.  The  Galerkin  finite 
element  approximation  is  then  gi'  n  by 

{findf/h  €  Tfcsuch  that  for  ea  V/(l  <  /  <  n) 

/n.  -  /,  ^U^)n^NId^H. 

(3.3) 

The  integrals  appearing  h  are  evaluated  in  the 
standard  finite  element  fc  by  summing  individual 
element  and  boundary  suriai.e  contributions,  the  com¬ 
pact  support  of  the  shape  function  Ni  means  that  the 
equation  can  be  written  as 


■NidQH  = 


y  ( 

E/  Fi{Uk)^dQ^-y  f  Fi{UH)niNidTK, 


Here  p,  p,  and  e  denote  the  density,  pressure,  and  spe¬ 
cific  total  energy,  respectively,  and  u,-  is  the  veloc¬ 
ity  of  the  flow  in  the  coordinate  direction  n.  This 
set  of  equations  is  completed  by  the  addition  of  the 
equation-of-state 


P=  (7-  l)p(«  -  ).  (2-3) 

which  is  valid  for  perfect  gas,  where  y  is  the  ratio  of 
the  specific  heats. 

In  the  sequel,  we  assume  that  fl  is  the  flow  do¬ 
main,  r  its  boundary,  and  nj  the  unit  outward  normal 
vector  to  the  boimdary. 

3.  VARIATIONAL  FORMULATION  AND 
FINITE  ELEMENT  APPROXIMATION 

Let  T  be  a  trial  function  space  and  W  a  weighting 
function  space,  both  defined  to  consist  of  all  suitably 
smooth  functions.  An  equivalent  variational  formula¬ 
tion  of  (2.1)  is  given  by 

find  [/€T  such  that  VW  €  W 

/n  -  fn  +  ir  Wdr  =  0. 

(3.1) 

Assuming  flh  is  a  classical  triangulation  of  (2  with  the 
nodes  numbered  from  1  to  n  and  Th  the  boundary  of 
Qa  ,  we  approximate  the  trial  and  weighting  spaces  T 
and  W  by  their  subspaces  of  finite  dimension  Th  and 
Wa,  which  respectively  are  defined  by 

Ta  =  |tfA(x,0  I  U^ix,t)  =  yU,(t)Ni{x) 

Wa  =  |H^h(*)  I  Wh{x)  =  ^o/Ar/(x)|  (3.2) 


(3.4) 

where  the  summation  extends  over  those  elements  e 
and  boundary  surfaces  b  that  contain  node  J.  Insert¬ 
ing  the  assumed  form  for  Uh  in  Eq.(3.4),  the  left-hand 
side  integral  can  be  evaluated  exactly  to  give 


where  Af  denotes  the  finite  element  consistent  mass 
matrix.  For  steady  state  computations,  M  can  be 
replaced  by  the  lumped  (diagonal)  mass  matrix,  de¬ 
noted  by  Ml. 

4.  EDGE-BASED  UPWIND 
FINITE  ELEMENT  SCHEME 

It  is  shown  in  the  appendix  that  for  any  interior 
node,  Eq.  (3.4)  can  be  written  as 


jrj 

[ML-^]i  =  yc\j(Fi+Fi)  (4.1) 
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where  m/  is  the  number  of  edges  connected  to  the 
node  /,  and 


cii  =  -y 

e€ll 


Qe  dNr 
4  dxj 


(4.2) 


in  3D.  The  coefficient  Cij  denotes  the  weight  applied 
to  the  average  value  of  the  flux  on  the  edge  that  con¬ 
nects  nodes  I  and  J,  to  obtain  the  contribution  made 
by  the  edge  to  node  7,  whereas  Cji  denotes  the  weight 


2 


applied  to  the  same  quantity  to  obtain  the  contribu¬ 
tion  made  by  the  edge  to  node  J.  It  can  be  easily 
verified  that  these  weights  possess  the  properties 


mt 

=  Q  for  all  /  ,  (4.3) 

IJ 

=  -C^ji  for  all  /  and  J  .  (4.4) 

For  notational  convenience,  we  define  the  vector  Ctj 
by  the  expression 


Cij  —  {Cjj,Cij,Cjj)  ,  (4.5) 

and  let  Ljj  denote  the  modulus  and  ib/j  denote  a  unit 
vector  in  the  direction  of  C/j,  then  Eq.  (4.1)  can  be 
written  as 

jrr 

[Ml  MFr  +  Fj)  =  Yl 

u  tJ 

where 

Ft  =  {F},Fj,Fj).kij  (4.7) 

Fj  =  {F),F],F^).ku  .  (4.8) 

The  alternative  procedure  for  obtaining  the  dis¬ 
crete  form  of  the  equations  is  now  apparent.  While 
with  the  element-ba^  data  structure  information  is 
gathered  from  all  the  nodes  of  each  element,  oper¬ 
ated  on  the  element,  and  then  scattered  back  to  the 
nodes  of  the  element,  the  edge-based  algorithm  gath* 
era  information  from  ail  the  nodes  of  each  edge,  op¬ 
erates  it  on  the  edge,  and  then  scatters  it  back  to 
the  nodes  of  the  edge.  The  property  of  conservation 
in  the  numerical  scheme  is  guaranteed  by  the  asym¬ 
metry  of  edge  coefficients  as  expressed  in  Eq.  (4,3). 
This  edge-based  data  structure  not  only  improves 
the  efficiency  of  the  algorithm  [1],  but  also  enables 
a  straightforward  implementation  of  upwind  schemes 
in  the  context  of  finite  element  methods.  It  is  clear 
that  Eq.  (4.6)  is  nothing  but  a  classic  Galerkin  fi¬ 
nite  element  scheme,  which  is  equivalent  to  a  cen¬ 
tral  difference  type  scheme.  By  using  the  results  of 
Eq.  (4.4),  this  scheme  allows  for  the  appearance  of 
chequerboarding  modes  and  thus  suffers  from  numer¬ 
ical  instabilities  unless  some  type  of  numerical  dissi¬ 
pation  in  the  form  of  artificial  viscosity  is  introduced. 
A  stable  scheme  can  be  constructed,  for  example  us¬ 
ing  Roe’s  flux  difference  splitting  [10],  to  replace  the 
actual  flux  function  Fjj  in  Eq.  (4.6)  by  Roe’s  numer¬ 
ical  flux  formula  Tij : 


with 


I  AA  1=1  Ai  1 


/  1  \ 
u 

V 

% 

\ 


Au  - 

Aw  -  k,A?t 
Atw  —  kxAft 
\  tiAu  -I-  wAw  -1-  wAtw  —  qk^qk  / 


|(4.11) 


I  Afl»,5  1  =  1  A4,5  I  ( 


Ap±^Aqk 

2c2 


ii±kgC 
v±kfC 
w  ±  k,BC 
\h±qkC 


(4.12) 

where  q^  =  u^-l-w^-l-u;*.  A?*  =  Auk^+Avky+Awk,, 
and  qk  =  ukr  +  +  ^k, .  The  bar  designates  Roe- 

averaged  quantities,  which  are  defined  by 


«  =  («/  +  +  '/pjjpi) 

V  =  {vi  ^>|■VJy/pJ|pl)/{l  4-  y/pjjpt) 

tw  =  (u>/  -t-  V)j\/pj/pi)l{\  -(-  \fp7fpi) 

h  =  {hi  +  hj\/pj/pi)/{l  +  \/pj/pi) 
c-it  -l){h-Q.h*q^)  . 


Furthermore,  the  eigenvalues  of  A  are  Ai  =  ft  and 
A4,5  =  ft  ±  c  . 

In  order  to  prevent  entropy  violation,  an  entropy 
fix  is  imposed.  When  an  eigen>^ue  A  reduces  to  zero, 
a  smoothed  value,  |  A  |*,  is  defined  to  replace  |  A  |. 


if  I  A  |>  e 
if  I  A  |<  e 


(4.13) 


where  e  =  KmaxiXj  —  A/,0).  Ff  is  a  small  constant. 

It  can  be  shown  that  this  scheme  is  equivalent  to 
the  first  order  finite  volume  upwind  cell-vertex  scheme 
based  on  a  dual  mesh.  There  are  many  different  ways 
to  achieve  higher  order  accuracy.  In  the  present  study, 
a  scheme  of  higher  order  accuracy  is  achieved  by  using 
upwind-biased  interpolations  of  the  solution  U  via  the 
MUSCL  approach  [llj.  This  leads  to  the  flux  function 


Fu  =  Ft  +  FJ-\  A{Ut,  UJ)  I  {UJ  -  Ut)  (4.14) 
where 


Fu  =  Fi^Fj-\Aij\{JJ]-Ui)  (4.9) 


where 


I  A,J  I  {Uj  -  Ui)  =1  AA  I  +  I  AA  1  +  I  AFi  I 

(4.10) 


Ft  =  F{Ut),  FJ  =  F(i;7)  .  (4.15) 

The  upwind-biased  interpolations  for  Ut  and  UJ  are 
defined  by 

Ut  =  Ur  +  J[(l  -  k)AJ  -Kl  -I-  k){Uj  -  U,)]  (4.16) 
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U7  k)^1  +  (1  +  kWj  -  Ui)]  (4.17) 

where  the  forward  and  backward  difference  operators 
are  given  by 

AJ=Ui-  U,.i  =  2{VU)i  •  1"  -  {Uj  -  U,)  (4.18) 

Aj  =  Uj+i  -Uj=  2{VU)j  ■  l^-'  -  {Uj  -  Ui)  (4.19) 

where  l^'^  =xj  —xi  is  the  length  vector  of  this  edge. 

The  parameter  k  can  be  chosen  to  control  a  fam¬ 
ily  of  difference  schemes  in  the  interpolation.  On 
structured  meshes  it  is  easy  to  show  that  k  =  —  1 
yields  a  fully  upwind  scheme,  Jb  =  0  yields  semi- 
upwind  approximation  (Fromm’s  scheme),  and  k  =  I 
yields  central  differencing.  The  value  k  =  1/3  leads 
to  a  third-order-accurate  upwind-biased  scheme,  al¬ 
though  third-order-accuracy  is  strictly  correct  only 
for  one-dimensional  calculations.  Nevertheless,  k  = 
1/3  was  used  in  the  calculations  presented  herein. 
With  higher  order  spatial  accuracy,  spurious  oscilla¬ 
tions  in  the  vicinity  of  shock  waves  are  expected  to 
occur.  Some  form  of  limiting  is  usually  required  to 
eliminate  these  numerical  oscillations  of  the  solution 
and  to  provide  some  kind  of  monotonicity  property. 
The  flux  limiter  modifies  the  upwind-biased  interpo¬ 
lation  Ui  and  Uj  and  the  Eqs.(4.16)  and  (4.17)  are 
replaced,  respectively,  by 


=  Ui  +  ^((1  -  kai)AJ  -(-  (1  -f-  kai){Uj  -  Uj)] 

(4.20) 

UJ  =Uj-  ^((1  -  kaj)Aj  -f  (1  -t-  kaj){Uj  -  U,)] 

(4.21) 

where  a  is  the  flux  limiter.  The  Van  Albada  limiter 
employed  in  this  study  acts  in  a  continuously  differ¬ 
entiable  manner  and  is  defined  by 


On  the  inflow  and  outflow,  a  characteristic  anal¬ 
ysis  based  on  the  ID  Riemann  invariants  is  used  to 
correct  the  computed  values  of  the  flow  variables  at 
the  time  step  n-(-l.  In  the  sequel,  the  '  -  indicates  the 
known  linearized  variables,  i.e.,  quar  5s  at  the  time 
step  n,  the  index  c  designates  the  co:  ited  variables 
at  the  time  step  n  4-  1,  and  the  it  ★  represents 

the  modified  variables  at  the  tim  p  n  4-  1  after 
applying  the  boundary  conditions 

Note  the  eigenvalues  and  char  eristic  variables 
are 


A  = 


/  \ 
Vn 
Vn 

Vn-kC 
\  V„  -  C  / 


,W  = 


\ 


(5.2) 


where  t;„  and  v,  =  (  )  are  normal  and  tangen- 

\Vt2j 

tial  velocity  components.  The  number  of  boundary 
conditions  that  has  to  be  imposed  is  equal  to  the 
number  of  negative  eigenvalues.  For  supersonic  in¬ 
flow  (t;„  <  — c),  all  the  eigenvalues  are  negative,  and 
therefore  all  the  variables  have  to  be  imposed.  In  this 
case,  all  the  variables  are  simply  reset  to  freestream 
values.  For  subsonic  inflow  (— c  <  <  0),  four  eigen¬ 
values  are  negative,  and  one  is  positive.  wi,W2,W3, 
and  u)4  are  defined  by  the  freestream  values,  while  ws 
is  determined  from  the  computed  state.  The  following 
equations  are  then  obtained; 


tT,*  =  v,oo 

<t'n*  +  ^=Vnoo  +  ^ 

~Vn*  4"  ^  —  ~^nc  4- 


(5.3) 


s/  =  mai{0. 


2AT{Uj-U,)+€ 
{Ajy  +  {Uj-Ui)^  +  €^ 


(4.22) 


aj  =  ma£{0, 


2A+{Uj  -U,)  +  ( 

{Ajy  +  {Uj  -  Cf/)2  4-  c 


} 


(4.23) 


where  e  is  a  very  smrdl  number  to  prevent  division 
by  zero  in  smooth  regions  of  the  flow.  Three  options 
exist  concerning  the  choice  of  interpolation  variables: 
conservative  variables,  primitive  variables,  ruid  char¬ 
acteristic  variables.  Using  limiters  on  characteristic 
variables  seems  to  give  the  best  results.  However,  the 
primitive  variables  are  used  in  this  study  for  the  sake 
of  computational  efficiency. 


By  combining  these  equations,  we  get  the  unknown 
variable  3 


I  ^n*  =  ^noo  +  ( pc  ^  (5-4) 

I  tT,*  =  V,oo 

(  P*  =  5[Poo  +  Pc  +  MVnoo  -  Wne)]  • 

For  subsonic  outflow  (0  <  Un  <  c),  only  one  eigen¬ 
value  is  negative,  and  pressure  is  then  imposed  by  the 
freestream  value.  wi,W2,tV3,  and  wa  are  determined 
from  the  computed  values.  The  following  relations 
are  then  obtained: 


5.  BOUNDARY  CONDITIONS 

The  treatment  of  boundary  conditions  is  very 
important  for  rapid  convergence  to  steady-state  and 
serves  as  non-reflecting  boundary  conditions  for  un¬ 
steady  computations. 

On  the  solid  walls,  the  normal  velocity  vanishes 
u„  =  0  .  (5.1) 


'p*=Pc  +  (2=^) 

<  Vn*  =  fne  +  i^^pc  )  (5-5) 

Vf*  *  Vtoo 
P*  —  Poo  • 

For  supersonic  outflow  (c  <  v„),  all  the  eigenvalues 
are  negative,  and  therefore  all  the  information  comes 
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from'  the  domain.  In  this  case,  nothing  needs  to  be 
imposed  and  all  the  values  are  computed  values. 

6.  TEMPORAL  DISCRETIZATION 

Equation(4.7)  represents  the  time  evaluation  of 
the  unknown  vector  Ui{t)  at  node  I  of  the  grid.  As¬ 
suming  that  the  nodal  values  (/"  are  known  at  time 
tn,  the  solution  is  advanced  over  a  time  step  At, 
to  time  t„+i  by  an  explicit  multi-stage  Runge-Kutta 
time-stepping  scheme  given  by 

=  uy 


i;(p)  ^  p=  l,2....,m 


■jn+l  ^  fj(m) 

with  the  parameters  Qp  assigned  appropriate  values. 
The  scheme  is  second  order  accurate  in  time.  For 
steady  state  computations,  implicit  residual  smooth¬ 
ing  and  local  time-stepping  are  used  to  accelerate  con¬ 
vergence  to  steady  state.  The  residual  smoothing  al¬ 
lows  the  use  of  larger  CFL  numbers  than  the  one  dic¬ 
tated  by  the  stability  of  the  original  scheme.  This 
is  accomplished  by  averaging  implicitly  the  residual 
with  values  at  neighboring  grid  points.  These  implicit 
equations  are  solved  approximately  by  using  several 
Jacobi  iterations.  The  local  time-stepping  uses  sepa¬ 
rately  a  maximum  allowable  step  size  for  each  node 
according  to  the  local  stability  analysis. 

7.  NUMERICAL  RESULTS 

All  the  grids  used  here  were  generated  by  the  ad- 
vemcing  front  technique  [12].  All  computations  used  a 
three-stage  Runge-Kutta  time-stepping  scheme  with 
local  time  stepping  and  implicit  residual  smoothing. 
The  computations  were  started  with  uniform  flow  and 
advanced  with  a  CFL  number  of  4.  The  L2  norm  of 
density  residual  is  taken  as  a  criterion  to  test  the  con¬ 
vergence  history. 

7.1  Channel  with  a  circular  bump  on  the  lower  wall 

The  first  test  case  is  the  well  known  Ni’s  test  case. 
It  is  a  transonic  flow  in  a  channel  with  a  10%  thick  cir¬ 
cular  bump  on  the  bottom.  The  length  of  the  channel 
is  3,  its  height  1,  etnd  its  width  0.5.  The  inlet  Mach 
number  is  0.675.  This  is  a  3D  simulation  of  a  2D  flow. 
This  simple  test  case  is  chosen  to  assess  both  accu¬ 
racy  and  convergence  of  the  numerical  scheme  and  to 
validate  the  implementation  of  the  code.  The  mesh, 
which  contains  13,891  grid  points,  68,097  elements 
and  4,442  boundary  points,  is  depicted  in  figure  lb. 
The  convergence  history  is  shown  in  Fig.  la,  where 
a  monotone  convergence  to  computer  machine  zero 
is  observed.  Fig.lc  displays  the  computed  pressure 
contours  in  the  flow  field.  The  Mach  number  distri¬ 
bution  on  the  lower  wall,  shown  in  Fig. Id  indicates 
that  there  is  only  one  grid  point  within  the  shock 


structure;  this  demonstrates  the  sharp  shock  captur¬ 
ing  ability  of  Roe’s  approximate  Riemann  solver  for 
the  solution  of  steady  problems. 

7.2  Wing/pylon/finned-store  configuration 

The  second  test  case  is  conducted  for  a 
wing/pylon/finned-store  configuration  reported  in 
reference  [13].  The  configuration  consists  of  a  clipped 
delta  wing  with  a  45  degree  sweep  comprised  from  a 
constant  NACA64010  synunetric  airfoil  section.  The 
wing  has  a  root  chord  of  15  in.,  a  semi-span  of  13 
in.  and  a  taper  ratio  of  0.134.  The  pylon  is  located  at 
mid-span  station  and  has  a  cross-section  chsiracterized 
by  a  flat  plate  closed  at  the  leading  and  trailing  edges 
by  a  symmetrical  ogive  shape.  The  width  of  the  pylon 
is  0.294  in.  The  four  fins  on  the  store  are  defined  by 
a  constant  NACA0008  airfoil  section  with  a  leading 
edge  sweep  of  45  degrees  and  a  truncated  tip.  The 
mesh  used  in  the  computation  is  shown  in  Fig. 2a.  It 
contains  274,953  grid  points,  1,518,770  elements  and 
33,046  boundary  points.  The  flow  solutions  are  pre¬ 
sented  at  a  Mach  number  of  0.95  and  sm  angle  of 
attack  of  zero  degree.  Figures  2b  and  2c  show  the 
pressure  contours  on  the  upper  and  lower  wing  sur¬ 
face,  respectively.  The  computed  pressure  coefficient 
distributions  are  compared  with  experimental  data  at 
two  spanwise  stations  in  Fig. 2d.  The  comparison  with 
experimental  data  is  excellent  on  both  the  upper  and 
lower  surface  up  to  70  percent  chord.  As  expected 
from  the  Euler  solution,  the  computation  predicts  a 
shock  location  which  is  downstream  of  that  measured 
by  the  experiment  due  to  the  lack  of  viscous  effect. 

7.3  ONERA  M6  Wing  configuration 

The  third,  well  documented  case  is  the  transonic 
flow  over  the  ONERA  M6  wing  configuration.  The 
M6  wing  has  a  leading  edge  sweep  angle  of  30  degree, 
^tn  aspect  of  3.8,  and  a  taper  ratio  of  0.562.  The  airfoil 
section  of  the  wing  is  the  ONERA  ”D”  airfoil,  which 
is  a  10%  maximum  thickness-to-chord  ratio  conven¬ 
tional  section.  The  flow  solutions  are  presented  at 
a  Mach  number  of  0.84  and  an  angle  of  attack  of 
3.06.  The  grid  adaption  scheme  was  used  for  this 
test  case.  The  final  adapted  mesh  contains  133,206 
grid  points,  738,669  elements,  and  17,155  boundary 
points  after  two  levels  of  refinement.  The  refinement 
of  high  gradient  regions  such  as  shocks,  leading  edge 
and  wing  tip  is  well  captured.  The  final  adapted  up¬ 
per  and  lower  surface  meshes  are  shown  in  Fig. 3a. 
The  pressure  contours  on  the  upper  wing  surface  and 
lower  surface,  Eire  displayed  in  Fig. 3b,  respectively. 
The  upper  surface  contours  clearly  show  the  sharply 
captured  lambda-type  shock  structure  formed  by  the 
two  inboard  shock  waves,  which  merge  together  near 
80%  semispeui  to  form  the  single  strong  shock  wave 
in  the  outboard  region  of  the  wing.  The  computed 
pressure  coefficient  distributions  are  compared  with 
experimental  data  [14]  at  four  spanwise  stations  in 
Fig.2c.  The  results  obtained  compare  closely  with 
experimental  data,  except  at  the  root  stations,  due  to 
lack  of  viscous  effects. 

7.4  F-18  fighter  configuration 

The  final  case  is  a  complete  F-18  fighter  con- 
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figuration,  which  includes  the  wing,  horizontal  and 
vertical  tails,  and  flow-through  engine  ducts.  The 
mesh,  which  contains  93,642  grid  points,  505,087  el¬ 
ements  and  15,421  boundary  points  for  the  half-span 
airplame,  is  shown  in  Fig.4b.  The  computations  were 
performed  at  a  free  stream  of  Mach  number  of  0.9  and 
an  angle  of  attack  of  3  degrees.  The  convergence  his¬ 
tory  is  depicted  in  flgure  4a.  The  solution  is  converged 
to  engineering  accuracy  (a  decrease  of  a  four  order-of- 
magnitude  in  the  L2  norm  of  the  density  residual)  in 
1300  time-steps.  It  required  a  total  of  10  CPU  hours 
on  a  single  processor  Cray  2.  The  computed  pressure 
contours  on  the  surface  of  the  airplane  are  shown  in 
Fig. 3c.  Some  of  the  features  occurring  in  this  flow 
regime,  such  as  the  canopy  and  wing  shocks,  are  well 
captured. 

8.  CONCLUSIONS 

An  edge-based  upwind  flnite  element  scheme  has 
been  developed  for  the  solutions  of  the  compressible 
Euler  equations  on  unstructured  grids.  The  numerical 
scheme  has  been  tested  and  validated  on  some  well 
documented  configurations.  An  example  application 
is  presented  for  a  complete  F-18  fighter  configuration 
to  demonstrate  the  accuracy  and  robustness  of  the 
proposed  algorithm. 
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APPENDIX 

In  this  appendix,  the  evaluation  of  inviscid  fluxes 
using  an  edge-based  data  structure  is  derived. 

RHS{I)  =  -  f  F^^dQH+  '  FifXjNidTH 
-  f  (nnode-2)Fl!  —dfln 

/  {nnode-2)FlN,  -dClh,  (1) 

JriK  i 

where  nnode  is  the  number  of  nodes  an  element, 
and  the  last  two  integrals  are  artific  ly  added  in 
order  to  derive  the  desired  formula.  Jy  using  the 
Green’s  formula,  we  obtain 

/  FjNf^dQH  =  i  /  Ffn^NiNidTH.  (2) 

JriK  dxj  2  Jpk 
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Note  that 


«  VT  •  < 


Using  £q.  (2),  Eq  (1)  can  be  written  as 
RHS(I) 

=  -  I  iF^  +  {rinode  -  2)FI Ni)^dQh 

JrtK  oxj 


[  Fifij 


Nidr^ 


+  ^  /  (nnode  -  2)FjnjNiNidrH  (3) 

For  an  interior  point,  the  boundary  integrals  can  be 
dropped  and  the  right-hand-side  becomes 

RHS{I) 

=  -  /  (/’>  (nnode  -  2)F\Ni)^dnH 

JtlK 

=  -  E  /  -  2)F/Arr)|^dnk 


=  -  E  /  (  E  -  2)^ 

fit-'n-  7^1 


=  -  E  /  (  E  -  l)FiNi)^dQH 

7ti  ^ 


Then  the  following  expression  is  obtained 
RHS(I) 

nnode 


=  -E(E/  +  (6) 

^•1  e6/ 

J^l 

which  can  be  further  simplified  as 

m, 

RHS(I)  =  ^C{j{Fi+F>),  (7) 


ij 


where  m/  is  the  number  of  edges  connected  to  the 
node  I,  2tnd 


I 


(8) 


(4)  in  3D. 


Fig. 2d.  Comparisons  between  computed  and  experimental  surface  pressure 
coefficient  for  the  wing/pylon/store  configuration  ;\/oo  =  0.95. o  =  0® 
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Fig.3a  Upper  and  lower  surface  meshes  for  M6  wing; 
nelem=738,669,tipoinsl33,206,nbounsl7,155 


I'lKl 


Fig.3b  Computed  pressure  contours  on  the  upper 
and  lower  surfaces;  Mqo  =  0.84,  a  =  3.06” 


Fig.3c  Compparison  between  computed  and  experimental 
surface  pressure  coefficient  for  M6wing 


Fig.4a  Convergence  History  for  the  F-18  Fighter 


Fig.4b  Surface  Mesh  of  Triangles  for  the  F-18  Fighter 
nelem=50S,087,  npoin=93,642,  nboun=15,421 


Fig.4c  Computed  Pressure  Contours  on  the  F-18  Fighter 
at  Moo  =  0.9,0  =  3® 


APPENDIX  4:  SAMPLE  RUN 


